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CHAPTER  I 


INTRODUCTION 

The  realization  of  machinr.  recognition  of  pictorial  data  has  long 
been  a  challenging  goal,  but  has  seldom  been  attempted  with  anything 
more  complex  than  alphabetic  characters.  In  this  thesis  the  task  of 
recognition  of  three  dimensional  objects  from  their  optical  images  is 
considered  from  the  view  of  identifying  and  estimating  the  translation 
and  rotation  of  the  object  with  respect  to  a  given  reference  frame. 

The  approach  presented  here  makes  use  of  the  theory  of  two- 
dimensional  moment  invariants  for  planar  geometric  figures  developed 
by  Ming-kue  Hu  [1],  Complete  systems  of  moment  invariants  under  trans¬ 
lation,  similitude  and  orthogonal  transformations  are  derived.  By 
carefully  utilizing  these  properties,  a  sample  set  is  constructed  in 
which  each  sample  is  represented  by  a  vector  which  characterizes  the 
image  for  a  certain  orientation  of  some  object  from  the  given  group. 

A  pattern  recognition  technique  is  then  described  in  which  a  parametric 
representation  of  the  input  signal  is  employed.  The  decision  process 
using  typical  samples  partitions  the  space  into  regions  that  envelop 
the  chosen  samples  of  a  class.  A  simulation  program  based  on  the  above 
outline  is  successfully  developed  which  not  only  identifies  objects, 
but  also  determines  their  o"ientation  and  position  in  space. 
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In  this  research  three  different  objects,  a  F-4B  PHANTOM  II, 
a  MIRAGE  III  C,  and  a  MIG  21  aircraft,  were  considered.  The  models  of 
these  aircraft  of  scale  1  :  72  are  mathematically  represented  by  a 
"wire  model"  structure  (2]  by  considering  a  finite  number  of  points, 
called  nodes,  on  the  body  of  the  modal.  These  nodes  are  Inter¬ 
connected  to  approximate  the  curved  and  planar  surfaces  on  the  body  of 
the  aircraft  by  a  number  of  straight  lines.  A  simulation  program 
is  used  to  generate  the  computational  model  of  the  silhouette  of  the 
aircraft  for  any  position  and  orientation  in  space.  Certain  Invariant 
properties  of  these  Images,  as  mentioned  earlier,  are  used  for  the 
construction  of  the  sample  sets.  Different  amounts  of  noise  were 
added  to  the  points  on  the  boundary  of  the  simulated  image  to  learn 
about  the  performance  of  this  method  when  using  unfocused,  hary  or 
unclear  optical  images. 

Besides  Che  Identification  of  aircraft,  this  approach  can  also 
be  used  for  various  other  problems.  A  television  camera  onboard  a 
docking  space  craft  can  cake  a  picture  of  a  docking  target  and  thus 
after  estimating  translation  and  rotation,  the  docking  craft  can 
position  itself  for  automatic  docking.  In  addition,  this  approach  can 
help  In  developing  a  robot  eye  for  use  on  automatic  assembly  lines  In 
Industries. 

A  review  of  different  pattern  recognition  techniques  related  to 
this  research  appears  in  Chapter  II.  The  problem  formulation  along 
with  a  complete  set  of  Invariant  moments  is  presented  in  Chapter  III. 
In  Chapter  IV,  the  computational  results  including  the  effect  of  noise 
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on  parameter  estimation  are  given.  Conclusions,  a  summary  of  the 
results,  and  future  research  areas  are  discussed  in  Chapter  V. 
Finally,  Appendix  1  presents  documentation  of  computer  programs  at 
the  end  of  the  thesis. 
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CHAPTER  II 

SURVEY  OF  PREVIOUS  WORK 

2.1  Introduction 

Munson  [3]  describes  the  pattern  recognition  process  in  terms  of 
the  following  three  stages  where  each  is  considered  as  an  independent 
component : 

TRANSDUCER  PREPROCESSOR  -»•  CLASSIFIER. 

However,  considering  the  first  step,  "it  is  doubtful  if  recognition 
occurs  before  the  eyes  are  directed  toward  a  (known)  object,  since 
otherwise  we  would  not  bother  to  look  at  the  object"  (4],  Thus  it 
would  seem  reasonable  to  use  the  raw  data  in  the  form  of  the  image  as 
bulk  memory  and  allow  the  transducer  to  search  for  "regions  of  interest" 
15].  In  considering  the  three  stages,  the  literature  overwhelmingly 
concentrates  on  the  various  aspects  of  classification.  It  is  here 
that  a  substantial  objection  can  be  raised.  Is  not  the  more  sig« 
nif leant  part  of  the  problem  that  of  characterizing  the  world  by  a  set 
of  properties  that  provide  the  desired  discrimination?  In  fact 
Seltridge  [6]  defines  pattern  recognition  solely  in  terms  of  "the 
extraction  of  significant  features  from  a  background  of  Irrelevant 
detail." 

Kazmierczak  and  Stelnbuch  [7]  state  that  "the  human  visual  system 
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is  capable  of  selecting  features  or  criteria  from  a  pattern  where  the 
statement  of  the  description  would  be  independent  cf  registration, 
skew,  size,  contrast,  deformation,  or  other  noise  effects."  What  is 
needed  according  to  Duda  [8]  are  "rugged  features".  "A  rugged  feature 
is  one  whose  presence  is  not  changed,  and  whose  characteristics  are 
not  greatly  altered,  by  normal  variations  in  the  image  of  a  character 
in  a  given  category."  It  is  emphasized,  and  this  is  an  important 
point,  that  no  general  theory  exists  to  allow  us  to  choose  what  fea¬ 
tures  are  relevant  for  a  particular  problem.  With  these  comments  in 
mind,  the  object  of  the  rest  of  the  chapter  is  to  present  a  discussion 
of  a  few  computer  methods  and  algorithms  used  in  conjunction  with  image 
analysis. 

2.2  Contour  Tracing 

One  of  the  approaches  to  reduce  the  amount  of  data  in  a  picture 
involves  scanning  a  picture  and  tracing  a  contour  or  outline  of  the 
figure  and  then  basing  the  recognition  or  classification  decision  cn 
this  information  (9].  It  is  well  known  [10],  that  "contours  carry  a 
ignificant  fraction  of  the  information  required  for  recognition  of 
image  objects."  Examples  of  this  approach  applied  to  character  recog¬ 
nition  are  discussed  in  the  literature  [11]  -  [13]. 

Hemami,  McGhee  and  Gardner  [14]  in  their  paper  presented  an 
algorithm  which  uses  the  information  contained  in  the  boundary  of  the 
pattern  by  successively  reading  the  coordinates  of  the  boundary  and 
developing  a  nonlinear  regression  analysis  technique  for  simultaneous 
estimation  of  rotation  and  translation  of  the  image  objects.  One  of 
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the  advantages  of  using  a  contour  description  is  that  the  latter  is 
independent  of  shape,  translation,  site  and  rotation  [15]. 

2.3  Conic  Section  Approximation 

Here  we  consider  shape  description  in  terms  of  conic  sections. 
An  individual  pattern  is  defined  as  a  non-negative  function,  f,  on  the 
real  plane,  subject  to  certain  constraints  on  position,  size,  orien¬ 
tation,  etc.  In  a  given  frame  of  reference  any  conic  section  may  be 
uniquely  represented  as 

Q(x,y)  ■  ax2  +  2hxy  +  by2  +  2gx  +  2fy  +  c  “  0,  (2-1) 

2  2  2  2  2 

where  a  +4h  +b  +4g  +4f  +  c  ■  1  and  the  first  non  zero  element 
of  the  vector  (  a,  2h,  b,  2g,  2f,  c)  is  positive.  Such  a  vector  will 
be  referred  to  as  a  conic  vector. 

The  true  Euclidian  distance  from  a  point  (u,v)  to  the  nearest 
point  on  the  conic  Q(x,y)  •  0  is  a  troublesome  quantity  to  evaluate, 
and  we  use  instead  the  quantity  |Q(uv)| ,  which  vanishes  on  and  only  on 
the  conic  and,  loosely  speaking,  takes  larger  values  for  points  (u,v) 
further  from  the  conic.  The  weighted  squared  discrepancy  between  the 
pattern,  P,  and  the  conic  may  then  be  defined  as 

D  -  // P(u,v)  |Q(uv)!2  du  dv,  (2-2) 

and  the  "best"  conic  is  that  for  which  D  is  the  least.  This  problem 
may  be  formulated  as  an  eigen  value  problem  of  order  six,  the  conic 
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vector  corresponding  to  the  smallest  root  defining  the  best  conic  [16]. 

The  best  conic  approximation  has  been  used  by  Paton  [16]  as  a 
discriminator  in  chromosome  analysis.  An  example  of  a  best  conic 
approximation  to  a  certain  chromosome  pattern  is  shown  in  Fig.  2.1. 


Figure  2.1  Line  Drawing  of  a  Chromosome  Pattern  and 
its  Best  Conic  Approximation 


2.  it  Shape  Descriptors 

Intuitively,  it  is  preferable  to  describe  an  object  using  gross 
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properties  rather  than  local  or  neighborhood  descriptors.  Shape  de¬ 
scription  is  undoubtedly  one  of  the  most  important  aspects  of  pattern 
recognition.  It  is  desirable  in  many  applications  to  describe  the 
structure  of  an  object  independently  of  orientation,  translation  or 
even  some  types  of  distortion.  Interesting  discussions  regarding 
recognition  of  shape  and  their  computer  models  can  be  found  in  papers 
of  Blum  [ 17 ] ,  ( 18] . 

Based  on  some  shape  description  theories,  an  interesting  shape 
descriptor  referred  to  as  a  medial  axis  transformation  (MAT)  was 
developed  by  Blum  [18].  He  describes  the  generating  model  which  is 
used  to  define  MAT:  "  Consider  a  continuous  isotropic  plane  that  has 
the  following  properties  at  each  point:  1)  excitation  -  each  point 
can  have  a  value  of  0  or  1,  2)  propagation  -  each  excited  point  excites 
an  adjacent  point  with  a  delay  proportional  to  the  distance,  and  3) 
refractory  or  dead  time  -  once  fixed,  an  excited  point  is  not  affected 
by  a  second  firing  for  some  arbitrary  interval  of  time.  A  visual 
stimulus  from  which  the  contours  or  edges  have  been  extracted  impinges 
on  such  a  plane  at  some  fixed  time  and  excites  the  plane  at  those 
points.  This  excitation  spreads  uniformly  in  all  directions  but  in 
such  a  way  that  the  waves  generated  do  not  flow  through  each  point." 

The  MAT  is  then  defined  as  the  locus  of  the  corners  in  the  wave- 
front.  The  (propagating)  contours  have  been  likened  to  the  front  of 
a  grassfire  ignited  on  the  pattern  boundary  and  the  MAT  is  then  the 
locus  of  points  where  the  fire  is  extinguished.  Several  examples  of 
MATs  are  shown  in  Fig.  2.2.  Note  shat,  if  the  MAT  turns  out  to  be  a 
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straight  line,  then  the  shape  of  the  object  under  consideration  is 
symmetrical.  The  gross  properties  of  the  transformation  are  obviously 
related  to  the  macrcscopic  and  structural  properties  of  the  pattern. 
This  is  demonstrated  by  Blum  [18]  using  both  a  sketch-like  represen¬ 
tation  of  a  human  and  its  distorted  version  as  shown  in  Figure  2.3. 

The  basic  properties  of  the  MAT  -amain  unchanged. 

2.5  Moment  Transformations 

A  set  of  two-dimensional  moment  invariants  have  been  found  by 
Hu  [1].  Based  upon  these  moment  invariants,  a  pattern-recognition 
theory  has  been  formulated  which  considers  two  patterns  to  be  similar 
if  they  differ  at  most  in  the  following  respects: 

(A)  Location 

(B)  Size 

(C)  Orientation 

With  the  help  of  moment  transformations,  we  can  find  a  pattern  function 
as  a  number  m  ■  F(P)  associated  with  each  pattern,  P,  whicn  is 

(a)  invariant  under  (A)  to  (C) ,  i.e.,  it  patterns  Pj  and  P2 
are  similar  according  to  (A)  to  (C)  then  F(P^)  -  F(P^); 

(b)  characteristic  of  dissimilar  patterns,  i.e.,  if  patterns 
P^  and  arc  not  similar  according  to  (A)  to  (C)  then  F(P^)  4  F(P2); 

(c)  easy  to  compute. 

If  such  a  function  cannot  be  found,  then  one  might  try  to  use 
several  functions  Fj(P) ,  F2(P),  ...,  which  satisfy  (a)  and  (c) ,  though 
not  (b),  in  the  hope  that  for  any  given  pair  P^  and  P2  of  dissimilar 
patterns,  at  least  one  of  the  functions  F^P)  would  give  F^fP^J^F^ (Pj) . 
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The  latter  is  the  approach  we  propose  to  follow. 

2.5.1  Definition  of  Two  Dimensional  Moments 

Let  there  be  N  points  equally  distributed  along  the  boundary  of 
a  certain  pattern  as  shown  in  Figure  2.4. 


Figure  2.4  Discrete  Representation  of  the  Boundary  of  a  Pattern 

Let  the  coordinates  of  these  points  be  (x^.yp,  (*2^2^*  •••»  ^xN,yN^ 
Then  the  two  dimensional  (p+q)**1  order  moments  are  defined  as 


m 


PQ 


1  N 

N  l  xiPyiq  ’  p  +  q  "  O.1*2.  •••  • 


(2-3) 


i-1 


•vtSi. 
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Also  let  (x,y)  be  the  centroid  of  the  given  pattern.  Then, 


N 

N  •  x  -  l  x* 
i-1 


1  N 

M  l  xi  *  m10 

i-1 


(2-4) 


Similarly 


y  “  ">oi  • 


2.5.2  Central  Moments 

The  central  momenta  are  defined  as 

“pq-  H 

i-1 


(2-5) 


(2-6) 


where 


x  "  m10  '  y  ’  m01  •  (2~7) 

It  should  be  noted  that  «  u^q  -  0.  Now,  let  us  consider  the  effect 
of  translation  on  these  central  moments.  Let 


x^'  -  Xj  +  a 

(2-8) 

yi'  =  yi  +  p 

(2-9) 
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where  a,  6  are  constants.  Then, 


1  N 

®10'  "  N  .I.xi' 


i-1 


1  p  N _ •_  a 

"  N  1*1  + 


1-1 


“10 


+  a 


(2-10) 


or 


x '  -  x  +  a 


(2-11) 


Similarly 


y'  -  y  +  B 


(2-12) 


Using  this  new  centroid  of  the  pattern,  let  us  now  calculate  the  new 
central  moments  under  translation. 


><,'  -  i  J1<Xi,-*,)P  (yi’-y,)q 


N 

1  V 


i  l  (x^+a-x-a)^5  (y^B-y- 
N  i-i 


B)* 


pq 


(2-13) 


Thus  the  central  moments  are  invariant  under  translation. 
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From  Eq.  (2-6)  it  is  quite  easy  to  express  central  moments  in 
terms  of  ordinary  moments.  For  the  first  three  orders , 


u00  “  “00  "  1 


(2-14) 


u  ■  u  *  0 

01  10 


(2-15) 


u  ■  a  -  (m  ) 

20  20  10 


(2-16) 


u  ■  m  -  (m  V 
02  02  01 


(2-17) 


u  ■  m  -  3m  m  +  2  (m  ) ' 
30  30  20  10  V  10 ; 


(2-18) 


U21  "  “21  “  “20“01  •  2“llm10  +  2{B10)  B01 


(2-19) 


x  ■  m  -mm  -2mm  +  2  (m^„ )  m,  ^ 

12  12  02  10  11  01  or  10 


(2-20) 


u02  "  “03  "  2“llm01  +  2(m01): 


(2-21) 


From  here  on,  for  simplicity  of  description,  all  moments  referred  to 
are  central  moments  and  will  be  simply  expressed  aa: 


u 


pq 


l^ 

N 


N 

l  *iP  Viq  • 

i-1 


(2-22) 
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2.5.3  Similitude  Moment  Invariants 


Under  a  similitude  transformation,  l.e.,  change  of  size,  we 


have 


x' 

y* 


u 


y  o 

0  y 


,  y  ”  constant  (2-23) 


Let  us  now  calculate  the  new  central  moments  u^q  after  the  transfor¬ 
mation. 


“pq 


N 


'  -  ^  ?  X,  'P  V,  ' 

P<!  N  A  i  yi 


(2-24) 


i-1 


|  l  S  Hv  r’  V 

N  i-1 


■  vP+q  M  y4* 


R 


i-1 


or 


*  p+q 

u  -yu 

pq  pq 


(2-25,) 


Therefore  we  have  the  following  absolute  similitude  moment  invariants: 


u 


Y 


I 

pq 

p+q 


upq  * 


p  +  q  -  2,  3,  . . . 


(2-26) 
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Using  similitude  invariants  of  central  moments,  pattern  identification 
can  easily  be  accomplished  independently  of  translation  and  size. 

2.5.4  Orientation  Moment  Invariants 

Under  the  orthogonal  transformation  of  rotation, 

x'  -  x  cos  6  -  y  sin  6  ,  (2-27) 

y'  -  x  sin  6  +  y  cos  6  .  (2-28) 


It  can  be  shown  l I]  that  the  three  second  order  moments  3atisfy  the 
following  relations: 

2u11'  -  (u2q  -  u02)  sin20  +  2UJJC.OS26  ,  (2-31) 


u2o'  +  u02 ' 


u20  +  u02 


(2-32) 
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(u2o' 


"02'»J  +  "V*1 


"  <“20  -  u02)2  +  4“ll2  '  <2-’3> 


There  are  two  ways  of  using  Eqs.  (2-31),  (2-32)  and  (2-33)  to  accom¬ 
plish  pattern  identification  independently  of  orientation: 

(.')  The  method  of  principal  axes:  if  the  angle  6  is  deter¬ 
mined  fret'  the  equation  (2-31)  to  make  u^1  ■  0,  then  we  have. 


tan  2e  -  -  2u^/  (u2q  “  u02^‘ 


(2-34) 


The  x* ,  y'  axes  determined  by  any  particular  value  of  6  satisfying 
Sq.  (2-34)  are  called  principal  axes  of  the  pattern.  With  added  re¬ 
striction,  suen  as  u^’  >  u^'  «id  u,^*  >  0,  0  can  be  determined 
uniquely.  Ttoments  determined  with  respect  to  <*uch  e  pair  of  principal 
exes  ore  independent  of  orientation. 

\B)  The  wet hod  of  orthogonal  moment  invariants:  the  two 
r-L^tloce  Eqa.  (2-32)  and  (2-33)  are  invariant  under  rotation,  and 
thay  can  be  uaeo  directly  for  or if ntation- independent  pattern  identi¬ 
fication.  Let  these  two  invariant  relations  he  called  M,  and  H-> 

L  4. 

respectively.  The  discrimination  property  can  also  be  increased  by 
including  higher-order  moment  invariants.  For  third-order  moments, 
we  ccn  show  tnat  tba  following  four  expressions  are  invariant  under 
orthogonal  transformation. 

-  (ujQ  “  3u^2^  ■*"  (3u2j^  —  u03^ 


(2-35) 
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M4  “  (u30  +  u12)2  +  (U21  +  u03)2  (2*36) 

M5  “  ^u30  "  3u12^  ^u30  +  u12^  ^u30  +  U12^2  "  3<u21  +  UQ3)21 

+  (3u2i  -  u03>  <u21  +  u03),13(u30  +  u12)2  "  (u21  +  u03)2’) 

(2-37) 

M6  “  (u20  ”  u02^  *(u30  +  u12)2  "  <u21  +  u03)^ 

+  4Ull  (u30  +  un)  (u21  +  uQ3)  (2-38) 

Similarly,  higher-order  orthogonal  moment  Invariants  can  be 
derived.  In  fact,  it  has  been  found  that  thnre  exists  a  complete 
system  of  infinitely  many  such  invariants  {!]. 

It  is  interesting  to  note  that  in  the  above  methods,  because 
of  complete  orientation  independence  property,  different  patterns 
which  could  be  obtained  from  each  other  by  just  proper  rotation,  such 
as  '6'  or  '9'  can  not  be  distinguished.  If  the  given  pattern  is  of 
circular  or  n-fold  rotational  symmetry,  then  the  determination  of  6  by 
Eq.  (2-34)  breaks  down.  This  is  due  to  the  fact  that  both  numerator 
and  denominator  are  zero  for  such  patterns. 

2.6  Sample  Set  Construction  and  Linear  Separatibility 

Based  on  the  features  extracted  from  the  optical  image  or  the 
pattern,  one  can  form  a  sample  set  which  characterizes  classes  from  a 
small  number  of  their  members.  In  addition  to  this,  pattern  recognition 
techniques  must  perform  a  basic  function  of  recognizing  a  new  input 


stimulus  and  classify  it  as  a  member  of  one  of  several  classes. 

Machine  learning,  the  automatic  accomplishment  of  classification,  re¬ 
quires  partitioning  the  vector  space  into  regions  so  that  each  region 
should  contain  mostly  members  of  a  single  class.  The  regions  so  con¬ 
structed  characterize  the  classes.  The  block  diagram  shown  in  Figure 
2.5  illustrates  a  general  pattern  recognition  system  that  exhibits  the 
functions  discussed  above  [19].  The  parameter  extractor  is  used  to 

decision 


Figure  2.5  General  Pattern  Recognition  System 

represent  the  machine's  environment  as  a  vector  in  an  N-d  linens  ional 
vector  space.  Machine  learning  is  employed  to  determine,  from 
sample  inputs,  the  best  method  of  partitioning  the  space  into  differ¬ 
ent  decision  regions,  and  the  decision  device,  which  implements  the 
regions  designed  through  machine  learning,  evaluates  new  input  stimuli 
and  classifies  them  according  to  the  region  in  which  they  are  con¬ 
tained. 

Most  of  the  machines  used  are  linear,  employing  only  linear 


discriminants,  i.e.  ,  correlations  with  stored  reference  vectors  or 
comparisons  of  weighted  combinations  of  the  parameters.  Hyper¬ 
planes  are  used  to  partition  the  space  of  measurable  input  parameters 
to  separate  members  of  one  stimulus  class  from  those  of  another. 
Combinations  of  these  linear  techniques  with  logical  rules,  can  con¬ 
struct  boundaries  to  quite  complex  distributions.  In  many  practical 
problems  the  classes  are  not  linearly  separable  (not  separable  by 
hyperplanes)  and  their  members  are  not  contained  in  disjointed  simply- 
connected  regions  of  space  of  observable  parameters.  In  these  cases, 
better  decision  rules  than  those  provided  by  linear  discriminants 
should  be  used  to  minimize  the  probability  of  decision  errors  [19]. 
This  is  illustrated  in  Figure  2.6  where  members  of  different  classes 
are  contained  in  regions  labeled  A  and  B.  These  classes  are  linearly 
separable  in  Figure  2.6a  but  not  in  the  more  complex  distribution 
shown  in  Figure  2.6b. 


(a) 

Figure  2.6  Different  Distributions  of  Two  Classes 
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2.7  Recognition  of  Three-Dimensional  Objects 

The  schemes  developed  by  many  authors  for  the  recognition  of  a 

solid  object  from  its  optical  image  require  the  perspective  trans¬ 
formation  of  a  three-dimensional  field,  with  hidden  lines  removed. 

The  perspective  projection  is  used  to  fit  the  given  picture.  The 
computer  oriented  techniques  of  forming  perspective  projection  of  a 
given  object  are  presented  by  Weiss  [20],  Comba  [21]  and  Loutrel  [22]. 

Guzman  [23]  and  Winston  [24]  in  their  work  in  the  field  of 
recognition  of  three  dimensional  objects  developed  a  scheme  which  re¬ 
cognizes  the  objects  irrespective  of  their  translations  and  rotations. 
They  consider  a  three-dimensional  structure  composed  of  bricks,  wedges 
and  other  simple  objects.  It  is  assumed  in  their  work  that  a  pre¬ 
processing  of  some  sort  has  taken  place,  and  the  picture  to  be  analyzed 
is  available  in  a  symbolic  format  of  points,  lines  and  surfaces.  The 
recognition  scheme  identifies  the  object  from  its  picture  by  selecting 
a  combination  of  surfaces  and  relating  it  to  an  object. 

The  only  two  pattern  recognition  schemes  in  the  literature,  to 
the  author's  knowledge,  which  deal  with  the  estimation  of  three  trans¬ 
lations  and  three  rotations  associated  with  the  object,  are  developed 
by  Roberts  [25]  and  Advanl  [26].  Roberts  [25]  in  his  work  assumes 
that  the  objects  seen  could  be  constructed  out  of  some  familiar  parts, 
called  the  models.  The  procedure  starts  by  first  converting  the  pic¬ 
ture  into  a  line  diagram;  then  the  points  in  the  line  diagram  which  fit 
a  transformation  of  some  model  are  found.  This  model  has  a  set  of 
topological  equivalent  points.  Finally  the  mean-square  error 


minimization  technique  with  some  threshold  is  used  to  eliminate  models 
which  fit  the  picture  topologically,  but  do  not  fit  exactly  without 
being  deformed.  This  scheme  does  not  yield  the  depth  information  but 
relies  on  a  support  theorem  for  the  purpose  of  estimating  the  depth 
or  the  translation  along  the  optical  axis  of  the  camera.  This  theorem 
requires  the  object  in  the  scene  to  be  supported  by  Lhe  ground  or 
another  object  resting  on  the  ground.  This  is  one  of  the  disadvantages 
in  Roberts'  recognition  algorithm. 

Advanl  (26]  developed  an  algorithm  to  estimate  the  three  trans¬ 
lations  and  three  rotations  of  an  object  from  its  silhouette  by  the 
use  of  regression  analysis.  Advanl,  in  his  work,  synthesizes  the  sil¬ 
houette  for  r  certain  translation  and  rotation,  and  then  tries  to 
match  this  synthesized  silhouette  with  the  given  cilhouette  of  the 
object  with  unknown  translation  and  rotation.  Advanl' s  method  is 
fairly  accurate,  even  in  the  presence  of  large  amounts  of  noise,  but 
the  main  drawback  of  this  algorithm  is  that  it  takes  a  much  longer  time 
for  recognition  than  what  would  be  needed  in  many  practical  appli¬ 
cations.  A  new  technique  for  recognition  of  three  dimensional  objects 
is  developed  in  this  research  which  has  the  potential  of  reaching  real 


time  identification. 


CHAPTER  III 


PROBLEM  FORMULATION 


3.1  Introduction 

The  problem  which  concerns  us  may  be  stated  as  follows:  a 
digital  computer  receives  an  optical  image  of  a  three-dimensional 
object,  and  on  the  basis  of  this  information  it  has  to  identify  the 
object  and  estimate  its  position  and  orientation  in  space. 

It  is  assumed  here  that  it  is  possible  to  obtain  the  silhouette 
of  the  picture  through  some  type  of  preprocessing  of  the  optical  image 
[5].  On  this  silhouette  a  number  of  equally  spaced  data  points  are 
generated  and  used  for  the  calculation  of  moments.  Using  these 
moments,  as  mentioned  earlier,  the  computer  should  assign  the  image 
to  a  certain  object,  and  estimate  the  six  parameters  to  be  defined 
later  in  section  3. A. 

3.2  Transformation  of  the  Real  World 

The  first  assumption  here  is  that  the  picture  is  a  view  of  the 
real  world  recorded  by  a  camera  or  other  comparable  device  and  there¬ 
for''.  that  the  image  is  a  perspective  transformation  of  a  three- 
dimensional  field.  This  transformation  is  a  projection  of  each  point 
in  the  viewing  space,  coward  a  focal  point,  onto  a  plane.  The  trans¬ 
formation  depends  on  the  camera  used,  the  enlargement  printing  process. 
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and,  of  course,  the  coordinate  system  the  reel  world  is  referred  to. 
Let  us  fix  the  real  world  coordinates  X,  Y  ana  Z  by  assuming  that  the 
focal  plane  is  the  Y  -  -2f  plane,  that  the  focal  point  is  at  X  ■  0, 

Y  ■  -f,  Z  -  0,  and  the  optical  axis  is  colinear  with  the  Y  axis.  In 
order  that  the  picture  not  be  a  reflection,  ve  choose  the  focal  plane 
in  front  of  the  camera.  Thus  the  focal  plane  is  really  the  plane  of 
the  print,  not  of  the  negative.  Let  U  and  V  represent  the  coordinates 
on  the  focal  plane  of  the  projections  for  the  points  in  the  reel 
world.  This  arrangement  is  shown  in  Figure  3.1. 


Figure  3.1  Camera  Transformation 


Let  (Xit  Y^,  Z^)  be  a  point  in  the  viewing  space  and  (Uj,  Vj) 
be  its  projection  on  the  focal  plane  as  shown  in  Figure  3.1  This 
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transformation  is  shown  below: 


”*  ■  vh  '  x‘  •  <J-» 

v‘  '  y~T7  '  M  »-2) 

where,  f,  is  the  focal  length  of  the  camera. 

3.3  Mathematical  Representation  of  Three-Dimensional  Objects 
with  the  Wire  Frame  Structure 

The  mathematical  representation  of  the  three-dimensional  object 
is  achieved  by  selecting  a  number  of  points ,  depending  on  the  com¬ 
plexity  of  the  structure,  on  the  body  of  the  prototype  model  of  a 
given  scale.  These  points  are  referred  to  as  nodes.  These  nodes  are 
then  appropriately  Interconnected  to  approximate  the  planar  and  curved 
surfaces  on  the  body  of  the  object  by  a  number  of  straight  lines.  This 
results  in  a  so  called  "wire-frame"  structure  [2],  because  the  straight 
lines  could  be  considered  as  wires  connected  between  the  fixed  points 
(nodes)  which  form  the  frame.  The  connections  between  different  nodes 
can  be  mathematically  expressed  by  a  connection  matrix  in  which  each 
row  represents  a  certain  node,  and  therefore  the  connection  matrix  has 
the  same  number  of  rows  as  the  nodes  selected  to  construct  the  wire¬ 
frame  structure.  The  different  column  entries  for  each  row  in  the 
matrix  show  the  nodes  connected  to  the  one  represented  by  the  row.  The 
trailing  zeroes  in  a  row  show  that  no  more  nodes  are  connected  to  the 
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node  represented  by  the  row.  One  such  example  appears  in  Figure  3.2. 
It  is  to  be  noted  here  that  the  nc**  row  of  the  connection  matrix, 
represents  the  nc^  node  of  the  object,  shown  in  Figure  3.2. 


Connection  Matrix: 
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Figure  3,2  Wire  Frame  Structure  and  Connection  Matrix 
for  a  Winged  Parallelepiped 

3.4  Transformation  Matrix 

Let  us  fix  a  new  coordinate  system,  xyz,  to  the  center  of  gravity 
of  the  object  to  be  identified.  The  coordinates  of  the  nodes,  selected 
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on  the  body  of  the  prototype  of  the  object  for  the  wire  frame  structure 
representation,  are  measured  with  respect  to  this  xyz  coordinate 
system.  Let  the  H  nodes,  selected  for  a  certain  object  to  be  identi¬ 
fied,  have  the  coordinates  [(x^,  y^,  z^),  i  ■  1,  2,  . . . ,  M] .  At  first, 
let  the  three  axes  x,  y  and  z  of  the  coordinate  system  fixed  to  the 
object  be  colinear  with  the  respective  axes,  X,  Y  and  Z  of  the  real 
world  coordinates.  This  is  shown  by  the  transformation  given  below  in 
Eq.  (3-3),  which  relates  the  coordinates  x^ ,  y^  and  z^  of  a  node  to  the 
real  world  coordinates  X^,  Y^  and  Z^. 


*1 

r 

xi 

*i 

• 

n 

zi 

zi 

(3-3) 


The  relative  orientation  of  two  arbitrary  orthogonal  systems  can 
be  specified  by  a  set  of  no  fewer  than  three  angles,  usually  called 
Euler  angles.  Although  the  concept  of  the  Euler  angles  is  universally 
used  in  several  applications,  there  Js  no  agreement  on  the  definition 
of  the  Euler  angles.  There  are  a  large  number  of  possible  choices  for 
the  three  angles  required  to  define  an  Euler  set  of  angles  [27].  A 
particular  set  of  Euler  angles  is  selected  here  with  the  aim  of  obtain¬ 
ing  moment  invariance  with  one  of  the  angles  of  the  set.  The  trans¬ 
formations  for  the  set  of  Euler  angles  selected  are  shown  in  Figures 
3.3  through  3.6.  The  symbolic  representation  used  in  these  figures  to 
accomplish  coordinate  transformations  is  straightforward  and  is 
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discussed  in  [28].  the  three  Euler  angles  9,  *  and  $  are  called  the 
elevation,  atiouth  a^d  roll  angles  respectively. 


x'  ■  -z  ain  6  +  x  cos  0 
y*  »  7 

z'  *2  cos  9  +  x  sin  9 
Figure  3.3  Elevation  Angle  Transformation 
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Figure  3.4  A2imuth  Angle  Transformation 


Figure  3.6  Symbolic  Representation  of  the  Euler  Angle  Sequence 


The  transformation  matrix,  which  relates  the  coordinates  of  a 
node  with  respect  to  the  two  different  systems  defined  earlier,  is 
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as  shown  below: 
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cos  9 

0 
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, 

-sin0 

0 
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0 

0 

1 

0 
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COS$ 

*i 

for  i  -  1,  2 . M  (3-4) 


After  giving  the  above  orientation  to  the  xyz  system,  let  us  translate 
the  origin  of  this  system  to  a  point  (A,  B,  C)  with  respect  to  the  real 
world  coordinate  system.  The  modified  transformation  is  given  below: 


xi 

CO80  COSlJ) 

-cose  sini  cosi 
+  sind  a  ini 

cose  s ini  a ini 
+  sin0  cosi 

xi 

A 

Yi 

- 

8  ini' 

cos i  cosi 

-cosi  sini 

n 

+ 

B 

zi 

-sine  cosi 

sine  sini  cosi 
+  CO80  sini 

-sine  sini  sini 
+  cos0  cosi 

H 

. 

C 

for  i  -  1,  2,  . . . ,  M  (3-5) 

With  the  use  of  this  transformation  matrix,  each  node  of  the 
wire  frame  structure  can  be  represented  in  the  XYZ  coordinate  system, 
given  the  three  translational  and  the  three  rotational  parameters. 

The  projection  of  these  M  nodes  can  be  obtained  on  the  image  plane  and 
used  to  construct  the  computational  model  of  the  silhouette  [26].  This 
method  of  simulating  the  silhouette  for  a  given  set  of  six  parameters 
is  used  later  for  the  sample  set  construction. 


r 
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3.5  Moment  Invariance  with  Angle  of  Elevation 

It  will  be  shown  here  that  the  silhouettes  obtained  for  different 
values  of  the  elevation  angle,  with  the  other  parameters  remaining  the 
same,  are  similar  in  size  and  shape  but  differ  only  in  rotation.  There¬ 
fore  the  moment  invariant  functions  dcrivec  in  section  2.5.4  character¬ 
ize  this  silhouette  for  the  other  given  parameters. 

It  will  be  assumed  that  the  optical  axis  of  the  camera  is 
directed  to  pass  through  the  center  of  gravity  of  the  object  to  be 
identified,  and  thus  the  parameters  A  and  C,  defined  earlier,  will  take 
smaller  values  compared  to  the  distance  of  the  object  along  the  optical 
axis.  This  can  be  accomplished  by  using  some  kind  of  feedback  system 
which  will  change  the  direction  of  the  optical  axis  so  as  to  bring  the 
centroid  of  the  silhouette  to  the  origin  of  the  UV  coordinate  system. 

The  projection  of  a  certain  point  in  the  viewing  space  onto  the 
focal  plane  is  given  by  the  following  transformation 


f  +  Y. 


*i 


f  +  Y. 


(3-6) 


(3-7) 


where  (U^,  V^)  are  the  coordinates  of  the  projected  point  in  the  focal 
plane.  With  the  assumption  that  the  parameters  A  and  C  are  small  in 
Eq.  (3-5),  we  have  the  following  relations: 
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'ii  -  (oii +  (coo*  eos*^  +  (-coo*  sin*)^  +  B  (3-8) 

X^  ■■  cosO(coo*)x^  +  (-cos8  sin*  cos*  +  sin8  sin*)yi 

+  (cos8  sin*  sin*  +  sin8  co8*)zj  (3-9) 

2.  ■  -8in0(coo*)Xi  +  (sinS  sin*  cos*  +  cose  sin*)y^ 

+  (~sin6  sin*  sin*  +  cos0  cos*)z£  (3-10) 

It  is  interesting  to  note  here  that  Y^  is  not  a  function  of  the 
eleva'ion  angle  6  and  hence  remains  the  sane  for  all  values  of  9, 
given  the  parameters  *,  *  and  B.  Thus  the  factor  K  “  f / (f  +  Y^)  in 
Eqs.  (3-6)  and  (3-7)  remains  constant  with  6,  for  the  given  values  of 
*,  *  and  B.  Therefore  we  have 


K  Xi  , 

(3-11) 

K  Yt  , 

(3-12) 

where  K  is  a  constant.  It  can  further  be  shown  that  the  expression 
2  2 

(X^  +  )  is  also  not  a  function  of  the  elevation  angle  6,  given  the 

assumption  A  *  C  ■  0.  Using  Eqs.  (3-9)  and  (3-10),  we  have 

Xj2  +  Z2  ■  cos2*  x,2  +  (sin2*  cos2*  +  sin2*)yi2 

+  (cos2*  +  sin2*  sin2*)z^2 


(3-13) 
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Therefore,  using  this  result  along  with  Eqs.  (3-11)  and  (3-12),  we 
have 

Ut2  +  -  K2(Xt2  +  Zi)2 

-  f«>,  4,  B)  (3-14) 

This  shows  that  when  the  elevation  angle  of  a  certain  object  is 
changed,  the  locus  of  the  projection  of  each  node  of  the  object  onto 
the  focal  plane  is  a  circle  with  the  center  as  the  origin.  Thus  we 
can  in  effect  say  that  the  result  of  varying  the  elevation  angle  of  a 
certain  object,  on  its  silhouette,  is  just  its  rotation  in  the  focal 
plane,  with  no  change  in  size  and  shape.  Perspective  projections  with 
different  elevation  angles,  but  with  the  same  azimuth  angle,  roll 
angle  and  distance  along  optical  axis,  appear  in  Figure  3.7.  The 
moment  functions  M^,  M2,  Hj,  M^,  M5  and  derived  earlier  in  section 
2.5.4  are  therefore  invariant  with  the  elevation  angle  6. 

M1  "  u20  +  u02  <3'15> 
M2  ”  (u20  ”  u02^  +  ^ull^  (3-16) 

“  ^u30  ”  3u12^  +  (3u21  “  u03^  (3-17) 

M4  “  (u30  +  u12^2  (u?.l  +  u03^ 


(3-18) 
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2  2 

-  (u30  -  3u12)  (u3Q  +  u21)  [ (u30  +  u12)  -  3(u21  +  uQ3)  ] 

+  (3u21  -  u03)  (u21  +  uQ3>  •  [3(u30  +  u12)2  -  (u21  +  u03)2] 

(3-19) 


"  (u20  ‘  u02>  ((u30  +  u12)2  “  <u21  +  u03)2] 

+  4u11<u30  +  u12)  (u2i  +  u03) 


(3-20) 


3.6  Moment  Invariance  with  Distance  along  Optical  Axis 

It  la  clear  that  as  an  object  is  moved  along  the  optical  axis, 
the  first  order  effect  on  the  image  is  just  a  change  in  size.  The 
second  order  effect  is  that  a  few  small  portions  of  the  image  may 
appear  or  disappear  when  the  object  is  coved  along  the  optical  axis. 
This  second  order  effect  diminishes  as  the  distance  of  the  object  from 
the  camera  increases. 

The  radius  of  gyration,  r,  of  a  planar  pattern  is  defined  as 
follows ; 


r  “  (U2Q  +  u02)is  (3-21) 

The  radius  of  gyration  is  directly  proportional  to  the  size  of  the 
image  or  inversely  proportional  to  the  distance  of  the  object  along 
the  optical  axis.  The  size  of  the  image  here  is  defined  as  the  mini¬ 
mum  radius  of  a  circle  required  to  enclose  the  given  image  completely 


with  its  center  at  the  centroid  of  the  image.  Thus  the  product  RK  of 
the  radius  of  gyration  of  the  image  and  the  distance  along  the  optical 
axis  of  the  object  is  a  constant. 


RK  ■  (U20  +  u02^  '  B 


(3-22) 


Therefore  the  radius  of  gyration,  r,  should  be  used  to  normalize  the 
moment  functions  of  Eqs.  (3-15)  through  (3-20)  to  obtain  the  similitude 
moment  Invariance.  Using  the  results  derived  in  section  2.5.3  con¬ 
cerning  similitude  moment  invariants,  we  have  the  following  moment 
functions  invariant  with  the  elevation  angle  6  and  the  distance  B  along 
the  optical  axis: 


"T  l(u20  “  u02)2  +  4ullZl 
r 


(3-23) 


M3’  -  -L.  ((u30  -  3u12)2  +  (3u12  -  u03)2]  (3-24) 

r 

V  "  h  ((u30  +  U12>2  +  <u21  +  u03)Zl  (3"25) 

r 

M5’  -  l(«30  -  3u12)  (u30  +  u21)  j_(u30  +  u12)2  -  3(u21  +  uQ3)2j 

+  Ou21  -  uQ3)  (u21  +  uQ3)  ^3(u3q  +  u12)2  -  (u21  +  v03)  j] 


(3-26) 
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"  U02J  ((u30  +  u12)2  '  <u21  4  u03)2} 


4  4u11(u30  4  u12>  <u21  4  u03>) 


(3-27) 


3.7  Elevation  Bias  Error 

Let  us  consider  a  silhouette  of  an  aircraft  for  a  certain 
orientation  as  shown  in  Figure  3.8.  Let  this  pattern  have  a  principal 
axis  as  shown  in  the  figure.  This  principal  axis  is  the  same  as  the 
perspective  projection  of  the  x  axis  of  the  coordinate  system  fixed 
to  the  ceuter  of  gravity  of  the  object  or  the  aircraft  in  the  present 
case.  Therefore  the  inclination  of  this  principal  axlB  of  the  pattern 
with  the  U  axis  in  the  projection  plane  is  the  elevation  angle  8.  The 
positive  direction  of  8  is  shown  in  the  figure. 


Figure  3.8  Elevation  Bias  Error  for  a  Pattern 
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The  inclination  6m  of  tha  major  axis  of  the  ellipse  fitted  to 
the  pat.em  can  be  found  from  Eq.  (2-34),  as  shown  below 

tan  26b  -  -2un/(u2o  -  u02)  (3-28) 

with  added  restriction,  such  as  ujg*  >  Ugj'  and  u^g'  >  0,  0B  can  be 
determined  uniquely  from  Eq.  (3-28). 

The  elevation  bias  error,  6#,  for  a  pattern  as  seen  from  Figure 
3.8,  is  defined  as  the  difference  between  the  elevation  angle  0  and 
the  inclination  0B  of  the  major  axis  of  the  ellipse  fitted  to  the 
pattern. 


6e  -  9  -  0B  (3-29) 

It  is  interesting  to  note  here  again  that  as  the  elevation  angle  is 
changed,  the  silhouette  of  the  object  just  rotates  in  the  projection 
plane  and  thus  the  elevation  bias  error,  6e,  remains  constant  for 
certain  given  values  of  che  azimuth  angle  $,  the  roll  angle  <j>  and  the 
distance  B  along  the  opticcl  axis.  Since  the  first  order  effect  on  the 
silhouette,  when  moving  the  object  along  the  optical  axis  is  just  a 
chonge  in  size,  we  can  say  that  the  elevation  bias  error  is  nearly 
invariant  with  the  distance  B  along  the  optical  axis. 

3.8  Sample  Set  Construction 

In  sections  3.5  and  3.6,  we  have  derived  five  different  fun¬ 
ctions  of  moments  of  the  planar  pattern,  which  are  invariant  with  the 
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elevation  angle  9  and  the  distance  B  along  the  optical  axis  of  the 
object  to  be  identified.  Thus,  these  five  functions  of  moments  can 
be  used  to  provide  a  parametric  representation  of  the  silhouette  of 
the  object  for  given  values  of  the  azimuth  and  roll  angles  in  a  five- 
dimensional  vector  space.  The  moment  functions  M2,  M-j,  M^,  and  Mg 
of  Eqs.  (3-23)  through  (3-27)  are  of  different  orders.  Invariant 
moment  functions  of  the  same  order  should  be  used  for  the  parametric 
representation  of  the  silhouette  in  order  to  give  equal  weight  or 
the  same  importance  to  each  parameter  used  In  the  recognition  pro¬ 
cedure.  Therefore,  let 


P1  " 

m2' 

(3-30) 

°2  " 

<m//6 

(3-31) 

P3  ’ 

(m49  4  /6 

(3-32) 

P4  “ 

(m594/12 

(3-33) 

°5  ' 

(M6-)4/8 

(3-34) 

and  then,  the  vector  (5  *  (p^,  Pj.  P^»  P^»  P5)  represents  the  given  sil¬ 
houette  in  a  five  dimensional  vector  space. 

A  simulation  program  for  obtaining  the  silhouette  of  a  given 
object  for  any  translation  and  orientation  [26]  is  used  here  for  con¬ 
struction  of  the  sample  set.  The  possible  range  of  variations  of  the 
azimuth  angle  \p  and  the  roll  angle  $  are  discretized;  and  then  for  all 
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of  the  different  combinations  of  these  discretized  values  of  the  two 

angles  with  the  values  of  A  and  C  taken  as  zero,  silhouettes  are 

generated  from  which  the  parameters  o^,  P3*  P^*  the  constant 

RK,  and  the  elevation  bias  error  6e  are  computed.  Thus  a  sample  set 

is  constructed  for  a  certain  object  and  then  repeated  the  same  way 

for  the  other  objects  to  be  identified. 

3.9  Nearest  Neighbor  Rule  for  the  Identification 
and  Estimation  of  ths  Parameters 

Let  the  vector  represent  the  ic^  pattern  in  a  five  dimensional 
vector  space  from  the  sample  set  (o^,  n^) ,  .....  (pn,  nn)  where  the 
n^'s  take  values  in  the  9et  £l,  2,  ....  K  Each  is  considered 
to  be  the  index  of  the  category  to  which  the  i**1  individual  belongs, 
and  each  pj  is  the  outcome  o'  the  set  of  measurements  made  upon  that 
individual  pattern.  For  brevity,  we  can  say  "p^  belongs  to  when 
we  precisely  mean  that  the  ic^  individual,  upon  which  measurements  p^ 
have  been  made,  belongs  to  category  n^. 

A  new  pair  (p,  n)  is  given,  where  only  the  measureuent  p  is 
raade ,  and  it  is  desired  to  estimate  o  by  utilizing  the  information  con¬ 
tained  in  the  sample  set  of  correctly  classified  vectors.  We  shall 
call 

e  el°l*  V  j 

a  nearest  neighbor  [29]  to  p  if 

min  dfo^  p)  «■  d(p’ ,  p)  i*l,2,  ...n.  (3-35) 
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The  unknown  Input  t,  by  the  nearest  neighbor  rule,  belongs  to 
the  category  or  the  class  n'  of  its  nearest  neighbor  o'  defined  in 
Eq.  (3-35).  A  mistake  is  made  when  n*  4  n.  T1T6  'nearest' neighbor  rule 
described  here  utilizes  only  the  classification  of  the  nearest  neigh¬ 
bor.  The  remaining  n-1  classifications  are  ignored. 

The  computations  involved  in  this  approach  are  fairly  simple. 

An  alternative  approach  using  linear  separatibllity  of  the  sample  set 
into  different  classes  could  possibly  be  used  for  the  identification, 
but  then  the  computations  Involved  would  be  more  complex  and  difficult. 
3.10  Summary 

In  this  chapter  we  first  considered  the  mathematical  repre¬ 
sentation  of  the  object  to  be  identified  and  of  the  optical  system. 

Some  functions  of  moments  which  are  invariant  with  the  elevation  angle 
0  and  the  distance  B  along  the  optical  axis  were  used  to  characterize 
the  silhouette  for  the  given  values  of  the  azimuth  angle  and  the 
roll  angle  A  procedure  was  then  described  to  form  the  sample  set; 
then  using  this  sample  set  a  nearest  neighbor  rule  useful  for  the 
identification  and  the  estimation  of  the  parameters  of  the  unknown 


object  was  discussed. 


CHAPTER  IV 


_  COMPUTATIONAL  RESULTS 

4.1  Introduction 

Wire  Frame  structure  representation  of  three  different  aircraft, 
a  F-4U  Phantom  II,  a  Mirage  IIIC,  and  a  MIG  21  is  presented  in  this 
chapter.  The  moment  functions  which  are  invariant  with  the  elevation 
angle  and  the  distance  along  the  optical  axis,  derived  in  Chapter  III. 
are  then  analyzed .  Using  the  paiametric  i epresentation  of  the  sil¬ 
houettes  generated,  the  construction  of  the  sample  sat  for  the  three 
aircraft  considered  is  shown.  The  results  obtained  by  the  application 
of  the  nearest  neighhor  rula  described  in  Chapter  III  for  the  identi¬ 
fication  and  estimation  of  parameters  ere  presented  later  in  this 
chapter.  The  effect  of  noise  on  identification  and  estimation  is  also 
studied. 

4.2  Wire  Frame  Structure  Representation 

Prototype  models  of  scale  1:72  are  used  for  the  purpose  of 
obtaining  the  wire  frame  structure  representation  of  the  three  differ¬ 
ent  aircraft  considered.  A  number  of  points  called  nodes  are  selected 
on  the  body  of  the  prototype  *o  represent  the  three  dimensional  objects, 
or  the  aircraft,  in  the  present  case.  The  number  of  nodes  required  and 
their  positions  on  the  body  of  the  prototype  for  a  complete  and  fair 
representation  of  the  three  dimensional  object  depends  on  the  object’s 
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complexity.  It  is  obviou3  that  the  number  of  nodes  necessary  to 
represent  a  curved  surface  would  be  more  than  that  necessary  to  rep¬ 
resent  a  planar  surface.  The  computation  time  required  to  generate 
the  silhouette  for  a  certain  given  translation  and  rotation  increases 
more  than  linearly  with  the  number  of  nodes  selected. 

The  coordinates  of  the  79  nodes  selected  on  the  Phantom  F-AB 
aircraft,  with  respect  to  the  center  of  gravity  of  the  aircraft,  appear 
in  Table  I.  The  connection  matrix  for  these  nodes  selected  is  given 
in  Table  II. 


TABLE  1 

COORDINATES  OF  THE  NODES  FOR  THE  PHANTOM  AIRCRAFT 
(With  Respect  to  Its  Center  of  Gravity  and  in  Feet) 


Node 

X 

y 

Z 

1 

4  .46 

-1  .37 

5.21 

2 

4.46 

1.37 

5.21 

3 

26,04 

-2.03 

-0.92 

4 

26.04 

2.03 

-0.92 

5 

26.04 

-2.03 

1  ,33 

6 

26.04 

2.03 

1  .33 

7 

18.32 

-2.96 

-2.20 

8 

18.32 

2.96 

-2.20 

9 

18.32 

-3.75 

-2.01 

10 

18.32 

3.75 

-2,01 

11 

17.85 

-3.75 

2.06 

12 

17  .65 

3.75 

2.06 

13 

1  7  .61 

-2  .59 

2.49 

14 

17.61 

2.50 

2.49 

15 

10.28 

-2.44 

-2.49 

16 

10.28 

2  .44 

-2.49 

17 

10  .26 

-4.29 

-1  .67 

18 

10  .29 

4.29 

-1.67 

19 

13.21 

-4.53 

1.8? 

20 

13.21 

4  .53 

1  .63 

21 

3 .50 

-4.53 

1.83 

47 


TABLE  2 — Continued 


Node 

Nodes 

Connected 

To 

1  60 

58 

56 

54 

0 

0 

3 

6  1 

7  C 

71 

72 

b  9 

C 

r. 

n 

i  62 

70 

71 

72 

5  C 

C 

0 

0 

I  63 

64 

3 

4 

5 

6 

65 

0 

i  69 

63 

66 

3 

b 

0 

C 

0 

65 

63 

5 

6 

68 

c 

n 

3 

66 

64 

67 

7 

8 

3 

4 

0 

67 

6  6 

69 

0 

0 

0 

0 

0 

!  68 

65 

76 

1  3 

1  4 

0 

0 

69 

67 

70 

7 

f 

15 

16 

0 

70 

69 

2 

61 

62 

4  1 

42 

c 

71 

73 

7b 

61 

62 

5  1 

5  2 

c 

72 

70 

73 

61 

62 

0 

0 

r\ 

O 

73 

72 

71 

C 

0 

0 

c 

0 

74 

71 

47 

b  8 

43 

44 

c 

c 

7  5 

76 

47 

48 

43 

44 

53 

59 

76 

77 

78 

1 

2 

68 

75 

n 

77 

76 

79 

1 

2 

0 

0 

0 

7e 

76 

79 

1 

2 

0 

0 

0 

79 

77 

78 

1 

2 

0 

0 

0 

The  trailing  zeroes  in  the  rows  of  the  connection  matrix  of 
Table  II  show  that  no  more  nodes  are  connected  to  the  node  representing 
that  row.  For  example,  node  77  in  only  connected  to  the  nodes  76,  79, 

1  and  2 . 

Similar  measurements  were  made  for  nodes  selected  for  the  other 
two  aircraft.  The  body  of  the  Mirage  is  less  -onplex  in  shape  than 
those  of  the  Phantom  and  MIG  aircraft,  and  therefore  it  is  only  nec¬ 
essary  to  select  59  nodes  for  the  Mirage.  The  coordinates  of  the  nodes 
and  their  connection  matrix  for  the  Mirage  and  MIG  aircraft  appear  be¬ 


low  in  Tables  III  through  VI. 


TABLE  3 


COORDINATES  OF  THE  NODES  FOR  THE  MIRAGE  AIRCRAFT 
(Wich  Respect  to  Its  Center  of  Gravity  and  in  Feet) 


Node 

X 

y 

z 

1 

23 .4e 

0.15 

-0.02 

2 

23  .48 

-0.15 

-0.02 

3 

18  .62 

1  .56 

-0  .02 

A 

18  .62 

-1  .56 

-0  .02 

5 

1  2  .80 

1  .50 

-0.95 

6 

1  2  .80 

-1  .50 

-0.95 

I  7 

12  .80 

2.76 

-C  .29 

j  8 

12.80 

-2.76 

-0.29 

9 

12  .80 

2  .76 

1.12 

10 

12.80 

-2.76 

1.12 

11 

12  .80 

1  .50 

1.51 

12 

12  .80 

-1.50 

1.51 

13 

5  .30 

1  .14 

-1  .76 

1  14 

5  .30 

-1.14 

-1  .76 

15 

6.1  1 

2.73 

-0.68 

1  6 

6  .1  1 

-2  .73 

-0  .68 

17 

6.71 

3.39 

1.57 

18 

6  .71 

-3.39 

1  .57 

19 

-2.17 

0.90 

-1  .58 

20 

-2  .17 

-0.90 

-1.58 

21 

-2  .59 

1  .74 

-1  .04 

22 

-2.59 

-1  .74 

-1  .07 

23 

-4  .51 

2.52 

2.38 

24 

-4  .51 

-2.52 

2.38 

25 

-10  .42 

12.63 

3.25 

26 

-10.42 

-12 .63 

3.25 

27 

-4.51 

2.52 

2.74 

28 

-4.51 

-2.52 

2.74 

29 

-11.11 

0.33 

-1.34 

30 

-11.11 

-0  .33 

-1  .34 

31 

-10  .60 

1  .89 

-0.23 

32 

-10.60 

-1 .89 

-0.23 

33 

-13.18 

3.09 

2 . 38 

34 

-13.18 

-3.09 

2.38 

35 

-12.49 

12.63 

3.25 

36 

-12  .49 

-12  .63 

3 

37 

-16  .93 

1  .62 

1  ... 

38 

-16  .93 

-1.62 

1  .81 

39 

-1 6  .93 

0.39 

-0.71 

40 

-16  .93 

-0.39 

-0.71 

4  1 

-20  .17 

1  .8? 

0.79 

42 

-20  .17 

-1.83 

0.79 

TABLE  3 — Continued 


Node 

1  - 

X 

V 

z 

4  3 

-20 .1  7 

1  .23 

1  .8  1 

44 

-20 .1 7 

-1.23 

1.81 

<♦5 

29 .03 

0  .00 

-0  .02 

*♦  6 

23.48 

0  .00 

-0.47 

47 

2  3  .4  8 

0.00 

-0.17 

48 

18.62 

0.00 

-1.34 

49 

18.62 

o.oc 

1 .39 

5C 

14.66 

0.00 

-2.93 

51 

12.80 

o.oc 

1  .93 

52 

10.76 

0.00 

-2 . 72 

53 

5.30 

0.00 

-2.39 

54 

-1  .72 

o.oc 

-1.97 

55 

-7  .66 

0.00 

-2.69 

56 

-17.89 

0.00 

-7.19 

57 

-21  .07 

o.co 

-7.19 

58 

-18.52 

0.00 

-1.37 

59 

-21.46 

0.00 

-1.07 

TABLE  4 


CONNECTION  MATRIX  FOR  THE  MIRAGE  AIRCRAFT 


Nodes  Connected  To 


45 

46 

45 

46 

48 

49 

48 

49 

50 

52 

50 

52 

1C 
1 1 
12 


4< 

4< 


TABLE  4 — Continued 


Node 


Nodes  Connected  To 


1  3 

53 

5 

15 

19 

52 

0 

0 

14 

53 

6 

16 

20 

52 

0 

0 

15 

7 

13 

17 

21 

0 

0 

0 

16 

8 

14 

18 

22 

0 

0 

0 

17 

51 

9 

15 

18 

23 

25 

27 

1  8 

51 

10 

16 

17 

24 

26 

28 

19 

54 

1  3 

21 

29 

0 

0 

0 

20 

54 

14 

22 

30 

0 

0 

0 

21 

15 

19 

23 

31 

0 

0 

0 

22 

16 

20 

24 

32 

0 

0 

0 

23 

17 

21 

25 

33 

35 

0 

0 

24 

18 

22 

26 

34 

36 

0 

0 

25 

17 

23 

27 

35 

0 

0 

0 

26 

18 

24 

28 

36 

0 

0 

0 

27 

17 

25 

28 

33 

0 

0 

0 

28 

18 

26 

27 

34 

0 

0 

0 

29 

5  4 

55 

58 

19 

31 

39 

56 

30 

54 

55 

56 

20 

32 

40 

56 

31 

21 

29 

33 

37 

41 

0 

0 

32 

22 

30 

34 

38 

42 

0 

0 

33 

23 

31 

34 

35 

37 

27 

0 

34 

24 

32 

33 

36 

38 

28 

0 

35 

25 

33 

23 

0 

0 

0 

0 

36 

26 

34 

24 

0 

0 

0 

0 

3  7 

31 

33 

38 

43 

0 

0 

0 

38 

32 

34 

37 

44 

0 

0 

0 

39 

58 

59 

29 

41 

40 

0 

0 

40 

58 

59 

30 

42 

39 

0 

0 

41 

31 

39 

43 

42 

0 

0 

0 

42 

32 

40 

l  4 

41 

0 

0 

0 

43 

37 

41 

44 

0 

0 

0 

0 

44 

38 

42 

43 

0 

0 

0 

0 

45 

46 

47 

1 

2 

0 

c 

0 

46 

45 

48 

1 

2 

0 

0 

0 

47 

45 

49 

1 

2 

0 

0 

48 

46 

50 

7 

4 

0 

0 

0 

49 

47 

51 

3 

4 

11 

12 

0 

30 

48 

52 

5 

6 

0 

0 

0 

51 

49 

1 1 

12 

1  7 

18 

0 

0 

52 

50 

53 

5 

6 

13 

14 

0 

53 

52 

54 

1  3 

14 

0 

0 

3 

TABLE  5 

COORDINATES  OF  THE  NODES  FOR  THE  MIG  21  AIRCRAFT 
(With  Respect  to  Its  Center  of  Gravity  and  In  Feet) 


o  o  o 


TABLE  6 — Continued 


Node 

Nodes 

Connected  To 

64 

63 

65 

3  7 

38 

39 

40 

0 

65 

64 

39 

40 

53 

54 

0 

0 

66 

67 

68 

43 

44 

53 

64 

0 

0  7 

66 

68 

C 

0 

0 

0 

0 

68 

66 

67 

69 

43 

44 

0 

0 

69 

68 

74 

23 

24 

0 

c 

0 

70 

71 

73 

2  5 

26 

0 

0 

3 

71 

70 

72 

25 

26 

0 

0 

0 

12 

71 

73 

25 

26 

0 

0 

0 

73 

70 

12 

25 

26 

0 

0 

0 

74 

69 

75 

1  1 

1  2 

0 

0 

0 

75 

74 

lb 

5 

6 

0 

0 

0 

lb 

75 

77 

76 

0 

0 

0 

c 

11 

76 

78 

0 

0 

0 

0 

0 

78 

76 

77 

79 

0 

0 

0 

0 

79 

55 

78 

1 

2 

0 

0 

4.3  Generation  of  Data  Points 

It  is  desired  to  represent  the  given  silhouette  of  a  certain 
object  by  some  finite  number  of  points  called  data  points.  These  data 
points  are  used  for  the  calculation  of  moments  for  the  given  silhouette 
The  most  realistic  situation  about  the  existence  of  these  data  points 
is  perhaps  to  have  them  uniformly  distributed  along  the  silhouette  or 
the  boundary  of  the  optical  image.  The  subroutine  that  generated  these 
data  points  is  called  DATA. 

Let  UVLE  be  the  distance  between  any  two  adjacent  data  points  on 
the  silhouette  or  the  boundary.  The  number  of  data  points  generated 
for  a  given  silhouette  is  inversely  proportional  to  the  distance  UVLE; 
and  therefoie  as  the  distance  UVLE  is  decreased,  the  computation  time 
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required  for  the  moment  calculations  is  increased.  However,  by  de¬ 
c-easing  the  distance  UVLE,  one  can  obtain  a  more  faithful  discrete 
representation  of  the  given  continuous  silhouette.  A  compromise  for 
the  value  of  UVl.F.,  on  the  basis  of  the  computation  time  and  the 
accuracy  desired,  is  made. 

For  a  certain  given  silhouette,  tha  variation  in  the  values  of 
the  parameter  03  of  Eq .  (3-32)  with  the  distance  UVLE  is  shown  in 
Figure  4.1.  The  parameter  Oj  varies  considerably  for  values  of  UVLE 
greater  than  0.0003  inch  and  remains  nearly  constant  for  values  of 
UVLE  less  than  0.0003  inch.  This  shows  that  for  values  of  UVLE  for 
which  any  of  the  parameters  p^,  02 ,  o-j*  0^  or  remain  nearly  constant, 

we  have  a  comparatively  more  exact  and  fair  representation  of  the  given 
continuous  silhouette.  A  compromise,  on  the  basi3  of  the  computer  core 
available,  the  computation  time  and  the  accuracy  desired,  is  made  for 
UVLE  equal  to  0.0002  inch. 


0  00002  0.0004  0.001  0  002  UVLE 


1 

Fleure  4,1  Variation  of  the  Parameter  0-  with  UVLF. 

-A 

I 

I 
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4.4  Moment  Invariance  with  the  Angle  of  Elevation 

It  was  shown  in  Chapter  III  that  the  vector  p  ■  (o^,  P2>  o^, 

04,  O5).  where  the  parameters  Pp  P2>  Pp  and  05  are  functions  of 
moments  which  are  invariant  with  the  elevation  angle  9  and  the  distance 
B  along  the  optical  axis,  characterizes  a  silhouette  with  given  azimuth 
and  roll  angles.  The  invariance  of  one  of  these  five  parameters,  Pp 
with  the  elevation  angle  is  shown  in  Table  VII.  Similar  results  hold 
for  the  other  parameters  also. 


TABLE  7 


VARIATION  OF  THE  PARAMETER  0,  WITH  THE  ELEVATION  ANGLE 
AIRCRAFT:  MIRAGE,  B  -  7000  ft,,  *  -  108 ,  4>  *  20° 


9 

in  degrees 

A  =  C 
*0  ft. 

A  =  C 
-50  ft. 

A  -  C 
-100  ft. 

-90.0 

0.8158 

C  .  8  1  4  2 

0.8142 

-75  .0 

0.8158 

C  .8  142 

0.8142 

-60 .0 

0.8158 

0.614b 

0.8146 

-45.0 

0.8158 

0.8153 

0.8153 

-30.0 

0.8158 

0.8158 

0.8158 

-15  .0 

0.6158 

0.8166 

0.8166 

0.0 

0  .81 58 

0.6171 

0.8171 

15  .0 

0.8158 

0.8175 

0.8175 

30  .0 

0.8158 

0.8177 

0.8177 

45.0 

0.8158 

0 .81 7S 

0.8179 

60  .0 

0.6158 

0.8183 

0.8183 

75  .0 

0.8158 

0.8182 

0.8162 

90 .0 

C.8158 

0.8180 

0 . 5  1  8  G 

- J 

It  is  seen  from  Table  VII  that  the  value  of  the  parameter  does 
not  vary  with  the  elevation  angle  9  for  the  case  when  the  center  of 
gravity  o:  the  object  to  be  identified  lies  on  the  optical  axis.  For 
the  other  two  cases  where  the  center  of  gravity  is  off  the  optical  axis, 
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there  is  a  slight  variation  of  lcaa  than  1*  in  the  value  of  thu 
parameter  For  the  nresent  case,  it  is  reasonable  to  assume  the 

maximum  value  of  the  off-axis  distances  A  and  C  to  be  about  30  ft.  only, 
as  for  values  greater  than  this  the  aircraft,  with  our  present  consid¬ 
erations,  will  not  be  completely  in  view.  This  is  sh-'wn  in  Figure  4.2 
where  the  silhouette  is  generated  for  the  off-axis  d i. t-  :ices  A  and  C 
equal  to  25  ft.,  and  the  optical  axis  distance  B  equal  to  10,000  ft. 

This  maximum  value  for  the  off-axis  distances  A  and  C  will  obviously 
decrease  when  smaller  optical  axis  distances  are  considered. 

4.5  Moment  Invariance  with  Distance  Along  the  Optical  Axis 

It  was  stated  in  the  last  chapter  that  the  radius  of  gyration  of 
the  optical  image  or  the  silhouette  is  inversely  proportional  to  the 
distance  of  the  object  along  the  optical  axis,  and  hence  the  product  RK 
of  the  radius  of  gyration  and  the  distance  along  the  optical  axis  re¬ 
mains  constant.  This  is  shown  to  be  valid  from  the  computational  re¬ 
sults  shown  in  Table  VIII. 

The  fact  that  the  product  RK  remains  nearly  invariant  with  the 
distance  along  the  optical  axis  vis  used  in  Chapter  III  to  normalize 

some  of  the  moment  functions  to  obtain  their  invariance  in  relation  to 

the  distance  along  the  optical  axis.  The  variation  of  one  of  the 

parameters  o;  with  the  distance  B  is  shown  in  Table  IX. 

4.6  Sample  Set  Construction 

A  simulation  program  which  appears  in  Appendix  I  is  used  to  gen¬ 
erate  the  silhouette  of  an  object,  given  its  three  translational  and 
three  rotational  parameters.  The  possible  range  of  variation  of  the 
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TABLE  8 


VARIATION  OF  kK  -  (u20  +  v^)*5  •  B  WITH  DISTANCE  B 
AIRCRAFT:  PHANTOM,  *  -  5® ,  6  «  75° ,  0  =  10® 


B 

In  feet 

A  -  C 
-0  ft. 

A  -  C 
-10  ft. 

A  -  C 
=20  ft. 

i  5000  .0 

27  38  .8 

2  738 .4 

2  737  .9 

5500 .0 

2739.9 

2739.4 

2738 .7 

6000 .0 

2740.2 

2739.9 

2739.2 

6500  .0 

2739.2 

2740.0 

2739.9 

7000  .0 

2740.0 

2739.6 

2741.0 

7500  .0 

2739.1 

2740.2 

274C .3 

8000  .0 

2740.2 

2739.8 

2739.6 

6500  .0 

2740.6 

2740.3 

2740 .3 

9000 .0 

27  39.1 

2739.1 

2738 .2 

9500  .0 

2739.2 

2738.9 

2738 .7 

10000  .0 

2740.7 

2740.4 

274C .2 

10500  .0 

2739.6 

2739.4 

2739 .1 

11000  .0 

2740  .3 

2740.0 

2740  .2 

11500  .0 

2741  .1 

2742.7 

2742 .5 

12000  .0 

2737.9 

2737.7 

2737.5 

TABLE  9 


VARIATION 

AIRCRAF 

OF  THE  PARAMETER  p* 
T:  PHANTOM,  <fi  -  5® , 

WITH  THE 
♦  -  75°, 

DISTANCE  B 

6  -  10° 

B 

A  -  C 

A  -  C 

A  -  C 

in  feet 

-0  ft. 

-10  ft. 

-  20  ft. 

5000.0 

0.2076 

0  .2077 

0 .20  81 

5500  .0 

0.2086 

0.2081 

0  .2C77 

6000.0 

0.2069 

0 .2087 

0 . 2083 

6500  .0 

0.2065 

0 .2068 

0 .2071 

7000.0 

0.2086 

0 .2086 

0.2087 

7500  .0 

0 .2082 

0  .2074 

0 .2080 

8000.0 

0.2078 

0.2077 

0  .2063 

8500.0 

0.2073 

0 .2072 

0.2077 

9000  .0 

0.2082 

0 .2  079 

0  .2079 

9500.0 

0.2065 

0  •  ?.  0  6  ^ 

0  .2064 

10000  .0 

0 . 2C95 

C  .2095 

0 .2095 

10500  .0 

0 .2063 

0 .2063 

j  .  2  0  6  2 

1  1000  .0 

0.2075 

0 .2074 

0.2081 

1  1500.0 

0.2071 

0.2074 

0  .  2C73 

12000*0 

0 .2085 

00 

r> 

rvj 

O 

0  .2084 
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azimuth  angle  and  the  roll  angle  4  are  0  to  360  degrees.  But  because 
of  the  symmetry  of  the  aircraft  about  their  major  axis  or  the  x-axia  of 
the  coordinate  system  fixed  to  the  center  of  gravity,  the  silhouettes 
for  the  four  roll  angles  d,  180  ,  180  +d  and  360  will  be  almost 

similar.  Thus  we  will  expect  the  parameters  r>^,  Ot,  03,  o^  and  05  to 
have  nearly  the  same  values  for  all  the  four  angles  mentioned  above; 
and  therefore  the  same  vector  0  “  (0^,  03*  c3*  0$)  will  represent 

the  four  silhouettes  for  these  roll  angles,  with  the  same  azimuth  angles. 


This  is  shown  in  the  Figure  4.3. 


Figure  4.3  Representation  of  the  Four  Different  Roll  Angles 
having  similar  Silhouettes 
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The  same  considerations  apply  to  the  azimuth  angle  also.  Making 
use  of  these  symmetric  properties  of  the  roll  and  azimuth  angles,  we 
may  consider  the  variations  of  $  and  4)  from  0  to  90  degrees  only.  The 
range  of  variations  of  the  azimuth  and  roll  angles  are  discretized  at 
an  interval  of  5  degrees  each;  and  then  for  all  of  the  different  combin¬ 
ations  of  these  di3cretized  values  of  the  two  angles  with  A  and  C  taken 
as  zero,  silhouettes  are  generated  from  which  the  parameters  p^  , 

Pi,  the  constant  RK,  and  the  elevation  bias  error  0e  are 

computed.  Thus  the  whole  sample  set  consisting  of  361  pattern  vectors 
for  each  aircraft  is  constructed.  Tables  X  and  XI  give  the  sample  set 
for  the  Phantom  an*  Mirage  aircraft  respectively.  The  variations  for 
*  and  t  in  these  tables  appear  in  steps  of  15  degrees,  but  these  tables 
were  computed  in  steps  of  5  degree  variations. 

4.7  Identification  and  Estimation  of  the  Parameters 

The  Nearest  Neighbor  Rule  described  in  Chapter  III  is  used  for 
the  identification  and  estimation  of  the  three  translational  and  the 
three  rotational  parameters.  The  sample  set  (p2»  n2)  ... 

(on>  n  )  is  stored  in  the  computer;  and  then  for  an  unknown  pattern 
vector  o,  its  nearest  neighbor  defined  in  Eq.  (3-35)  is  searched.  This 
is  done  by  using  the  subroutine  ESTM,  given  in  Appendix  I. 

4.7.1  Identification  of  the  Unknown  Object 

Let  the  nearest  neighbor  of  the  unknown  pattern  vector  o  be 
(o',  n').  We  may  therefore  say  that  p  belongs  to  n'.  The  information 
contained  in  n'  is  the  class  of  the  object,  the  azimuth  and  roll  angles. 
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SAMPLE  SET  CONSTRUCTION  FOR  PHANTOM 
A  -  C  =  0,  B  =•  7000  feet 
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the  constant  RK  and  the  elevation  bias  error  6e.  Thus  the  class  of  the 
object,  or  the  aircraft  type  in  our  present  case,  obtained  from  n’  does 
the  required  identification.  We  define  a  new  term  certainty  of  identi¬ 
fication  as 


CER  -  D2/(D1  +  D2) 


(4-1) 


where 


“  distance  of  unknown  pattern  vector 
from  its  nearest  neighbor  of  the 
class  identified 

D2  ■  distance  of  unknown  pattern  vector 
from  its  nearest  neighbor  of  the 
ether  class. 

4.7.2  Estimation  of  the  Translational  and 
Rotational  Parameters 

Information  about  the  azimuth  and  roil  angles,  the  constant  RK 
and  the  elevation  bias  error  0e  is  obtained  from  n' ,  and  is  used  to 
estimate  the  six  parameters. 

The  values  of  the  azimuth  and  roll  angles  obtained  from  n’  are 
only  approximate  estimates  for  the  two  angles.  A  technique  using  linear 
interpolation  is  then  used  to  generate  102  new  pattern  vectors  of  some 
known  azimuth  and  roll  angles  around  the  nearest  neighbor  (o’,  n'). 

Let  dig  and  pg  be  the  azimuth  and  roll  angles  obtained  from  n'.  The 
closest  8  neighbors  of  o'  from  the  sample  set  are  the  pattern  vectors 
with  the  azimuth  and  roll  angles  as  shown  below  in  Figure  4.4.  Using 
''hese  8  neighbors  of  o',  we  can  Generate  102  new  pattern  vectors  for 


(tc 


Figure  4. A  Values  of  ^  arid  $  for  New  Pattern 
Vectors  Generated  by  Interpolations 
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the  variation  in  steps  of  1  degree  each  in  the  values  cc  the  azimuth 

and  roll  angles.  This  is  accomplished  by  interpolating  the  values  of 

o.  ,  c0,  o,.  o,  and  c,  senarate'v  for  all  the  new  positions  shown  in 
12  3«  5' 

Figure  4.3  by  using  Che  above  mentioned  8  neiehbors.  This  is  done  by 
the  INTP  subroutine.  The  Nearest  Neighbor  Rule  is  then  used  again  to 
find  the  closest  match  for  the  unknown  vector  <$  from  the  newlv  gener¬ 
ated  pattern  vectors.  Thus  from  this  closest  match  found,  we  extract 
estimates  for  the  azimuth  and  roll  angles  accurate  up  to  1  degree  in 
noiseless  cases. 

The  estimate  for  the  distance  of  the  unknown  object  along  the 
optical  axis  is  as  shown  below: 


Be  -  RKjr/r  (4-2) 

where  R«E  is  the  value  of  the  constant  RK  obtained  from  n' ,  and  r  is 
the  radius  of  gyration  of  the  silhouette  of  the  unknown  object  to  be 
identif led. 

The  inclination  of  the  major  axis  of  the  ellipse  fitted  to 
the  pattern  can  be  found  from  Eq.  (3-28),  which  appears  below 
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The  elevation  bias  error  6fi  for  the  unknown  pattern  vector  ia  obtained 
from  n' ,  and  therefore  the  estimate  6g  for  the  elevation  angle  is 
given  by 

eE  "  em  +  ®e  (A--4) 

For  the  purpose  of  estimating  the  values  of  the  off-axis 
distances  A  and  Ct  we  assume  chat  the  center  of  gravity  of  the  unknown 
object  is  projected  very  near  the  centroid  of  its  silhouette  in  the 
image  plane.  Therefore,  from  Eqs.  (3-1)  and  (3-2),  we  have, 

aE  *  XCG  *  UCC  '  (Be  +  f>/£  (4-5) 

CE  “  YCG  "  VCG  ‘  <BE  +  f>/f  <4"6> 

where  (uCq>  *s  the  centro^  of  t*ie  silhouette,  B£  is  the  estimate 

for  the  distance  along  the  optical  axis  and  f  is  the  focal  length  of 
the  optical  system.  Thus  we  have  been  able  to  estimate  the  three  trans¬ 
lational  and  the  three  rotational  parameters  of  the  unknown  object. 

4.7.3  Results  of  Identification  and  Estimation 
in  the  Absence  of  Noise 

Several  silhouettes  were  generated  with  known  reference  trans¬ 
lational  and  rotational  parameters,  and  the  pattern  vector  p  ■  (p^,  p j, 
P3,  o^,  D5)  was  computed  for  each  of  the  silhouettes.  With  the  help  of 
these  pattern  vectors,  the  identification  and  estimation  was  achieved. 
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Results  of  some  of  the  identifications  and  estimitions  appear  in  Table 
XII. 


It  is  seen  from  the  results  of  Table  XII,  that  the  accuracv  of 
the  estimation  of  parameters  is  as  shown  below  in  Table  XIII. 


TABLE  13 

ACCURACY  OF  ESTIMATION  OF  THF.  PARAMETERS 


Parameters 


Accuracy 


Azimuth  Angle 
Roll  Angle 
Elevation  Angle 
Distance  B 
Distance  A 
Distance  C 


±1  degree 
±1  degree 
♦0.1  degree 
2  percent 
±2.0  feet 
±2.0  feet 


It  is  to  be  noted  from  Table  XII  that  the  certainty  of  Identi¬ 
fication  was  always  found  to  be  very  high,  and  from  this  we  can  Infer 
that  the  clusters  formed  by  the  pattern  vectors  of  each  aircraft  are 
nonover lapping  or  separate. 

4.7.4  Results  of  Identif icaticn  and  Estimation 
in  the  Presence  of  Noise 

In  a  realistic  system  one  would  expect  that  the  measurement 
points  would  not  exactlv  overlay  the  boundary  of  the  pattern  from 
which  they  came.  This  error  may  be  due  to  several  different  reasons. 
For  instance,  if  the  field  of  view  has  been  slightly  clouded  over,  or 
defocused,  then  the  boundary  of  the  pattern  is  no  longer  precise  and 
the  exact  coordinates  of  points  lying  on  the  boundary  can  onlv  he  esti¬ 
mated.  Even  if  the  field  of  view  is  clear  there  is  still  the 
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possibility  that  the  electronic  equipment  associated  with  the  optical 
system  can  commit  errors,  he  they  internal  or  transmission  errors. 
Furthermore,  there  is  always  the  quantization  error  associated  with 
analog  to  digital  conversion.  All  of  these  errors  may  be  considered 
as  forms  of  noise. 

A  simulated  silhouette  for  certain  reference  parameters  for  the 
Mirage  appears  in  Figure  4.5  with  a  certain  amount  ot  noise  added  to 
the  data  points.  The  identification  and  estimation  for  this  simulated 
silhouette  is  shown  below  in  Table  XIV. 


TABLE  14 

IDENTIFICATION  AND  ESTIMATION  IN  THE  PRESENCE  OF  NOISE 
(Aircraft  Considered:  Mirage,  Aircraft  was  Identified  to  be  Mirage) 
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Simulated  silhouettes  for  certain  reference  parameters  for  the 
Phantom  aircraft  with  different  amounts  of  noise  added  appear  in  Figure 
4.6  through  4.10.  The  identif ication  and  estimation  for  all  these  sim¬ 
ulated  silhouettes  appear  in  Table  15.  It  should  again  be  noted  that 
the  identification  of  the  aircraft  type  was  done  correctly  even  in  the 
presence  of  noise,  and  therefore  showing  the  presence  of  nonover fanning 
clusters  formed  by  the  pattern  vectors  of  each  aircraft. 


7  Perspective  View  of  Phantom  wi 
ft.,  C  -  12  ft.,  <1/  «  7*.  *  »  22’ 
!lse  of  2.2%  added  to  the  data  po: 
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TABLE  15 

IDENTIFICATION  AND  ESTIMATION  IN  THE  PRESENCE  OF  NOISE 
AIRCRAFT  CONSIDERED:  PHANTOM 

REFERENCE  PARAMETERS:  A  -  10  FEET,  B  -  7500  FEET,  C  -  12  FEET, 
■i  ■  7*.  *  -  22*.  6  -  18° 

(Aircraft  was  always  identified  to  be  Phantom) 


Noise 

A 

B 

Estimated 

C 

Parameters 

6 

CER 

0.0% 

3.6 

7468 

12.3 

8 

22 

18.1 

77.8 

2.2% 

8.5 

7459 

22.0 

7 

22 

18.1 

82.8 

4.8% 

8.6 

7512 

12.4 

5 

23 

18.3 

86.4 

7.6  % 

8.6 

7474 

12.4 

7 

22 

18.1 

81.9 

10.0% 

8.7 

7479 

12.3 

2 

24 

18.2 

94.6 

CHAPTER  V 


SUMMARY  AND  CONCLUSIONS 

This  thesis  is  aimed  at  a  solution  to  the  problem  of  real  time 
identification  of  three  dimensional  objects  from  their  optical  Images. 
The  approach  taken  here  relies  upon  the  invariant  functions  of  moments 
which  are  used  for  the  construction  of  the  sample  set. 

The  method  of  moments  appears  to  work  quite  well  for  the  identi¬ 
fication  of  aircraft  from  their  optical  images.  The  results  obtained 
in  Chapter  IV  show  that  the  performance  of  the  method  of  moments  for 
aircraft  Identification  is  comparable  to  that  of  a  human  photo  inter¬ 
preter,  and  it  is  roughly  equivalent  to  or  better  than  heuristic  tech¬ 
niques.  Failure  to  surpass  the  performance  o  a  human  is  not  always  a 
serious  liability,  however.  In  many  situations,  automation  is  re¬ 
quired,  even  if  the  resulting  performance  is  not  as  good  as  that  of  a 
human.  When  automation  is  required,  the  method  of  moments  is  quite 
possibly  the  best  method,  because  it  is  easier  to  implement  than  the 
other  heuristic  techniques,  which  require  point-by-point  image  analysis. 

So  far  as  is  known  to  the  author,  no  alternative  Image  pro¬ 
cessing  technique  exists  with  a  capability  of  real  time  identification 
and  estimation  of  parameters  of  three  dimensional  objects.  The 
algorithm  developed  in  this  research  has  the  potential  of  reaching  real 
time  identification  when  a  special  purpose  computer  with  multiprocessing 
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capabilities  is  used  to  do  the  requisite  searches.  An  associative 
memory  could  also  be  used  to  store  the  sample  set  in  order  to  de¬ 
crease  the  identification  time  required. 

While  all  of  the  results  contained  in  this  thesis  relate  to 
three  different  aircraft  that  were  considered,  it  appears  that  this 
approach  would  be  applicable  for  a  larger  number  of  objects  as  well. 
Further  work  in  the  following  areas  could  lead  to  better  accuracy  in 
estimation  of  the  parameters  and  in  reduction  of  the  limitations  of 
the  present  method: 

(1)  An  optimum  distance  function  for  use  in  the  Nearest  Neighbor 
Rule  described  in  Chapter  III. 

(2)  Linear  separatlbility  for  dividing  the  sample  set  into  dif¬ 
ferent  regions  or  classes. 

(3)  Syntactic  analysis  to  resolve  the  ambiguities  between  certain 
roll  angles  and  azimuth  angles  such  as  between  if  and  -if,  if  and  -if  or  the 
nose  and  the  tail  of  the  aircraft. 

(4)  Design  of  a  system  to  gather  data,  for  use  in  moment  cal¬ 
culations,  from  the  optical  image  of  an  object. 


APPENDIX  I:  COMPUTER  PROGRAMS 


o  i M  e  s  s  i  a  n 

CL!6).CM6).MC0N(80,7)  ,U0(  3000)  ,VD<  3000)  , 

1  C  1  I  6) ,C  <6 ) .P  <80) ,0!  80) ,R( 80 ) 

DIMENSION  U8N0(80 ) ,V6ND! 60) 

DIMENSION  RM6  (  1  .  1  )  «RM  1  .  1  )  . 

1  R*  1  (  1  . 1)  ,RM2(  I  .1 )  ,RM3(  1  , 1  )  ,RM4<  1  ,  1  >  ,RM5(  1  .  1  )  , 
1  E6(  1  .  1  I 

RE AQ(6 . 100 1 )  F.RPAT.UVLE.NPAR.NOP 

1001  FORMAT!  F6.1.F10.2.F5.1.2I3) 

REACH  6  1002  )  (P(  I  )  .Q  (  I  )  » R  (  N  #  I  *  1  «NQP  ) 

1002  FORM4T ( 3F9.4  ) 

REA0(6.1003)  CO (1)  ,  CO < 2 )  , C 0 t 3  ) 

1003  FORMAT ( 3F7 .1 ) 

READ! 6. 1006)  ( (MCOM l , 1 1 ) , 1 1* 1 .7) . I =1 ,NQP  ) 

1006  F0RMAK713) 

CALL  RELEAS  (6) 

F  =  12  . 

UV  LE  =  G  .0002 
TYPE  6C0 

600  FORMAT!  '  PLEASE  TYPE  OEC  TAPE  NO.  IN  FORMAT  11 
RE  AO (5.8886)  IU 
8886  FORMAT ( 1 1 ) 
l  U  *  1  U  ♦  8 
TYPE  2211 

2211  FORMAT!  '  ENTER  AIRCRAFT  CODE  •) 

REA0I5.2213)  PLANE 
2213  FORMAT ( F2.0  ) 

TYPE  bO  3 

603  FORMAT!  •  ENTER  CO-ORD  INATf.S  X,ZY  IN  FORMAT 

2 ( F5. 1  . 1  X  ) 

1  . F8 .1  '  ) 

read < 5 . 604 )  coil) , com, com 

604  FORMAT !2(F5 .1 .1X1  »F8  .1) 

TYPE  601 

bOl  FORMAT!'  ENTER  ROLL.ASYMTH  ELEVATION  ANGLES  !N 
FORMAT 

13! F5  .1  .IX)  '  ) 

REAJ(5  .602  )  CO! 4 ) ,C0 ( 5  )  ,C0 ( 6  ) 

602  FORMAT ( 3! F5 . 1  .  IX  )  ) 

TYPE  bO  5 

6C  5  FORMAT!'  ENTER  SO  IN  FORMAT  Fb.4  ') 

R  E  AO ( 5  .606  )  SO 
6Cb  FORMAT  ( F6 .4  ) 

TYPE  1245 
RE  AD!  5 . 2468  )  1FN 
2468  FORMAT (A3) 

1  245  FORMAT!'  ENTER  FILE  NAME  AS  ',/) 

CALL  OF  I LE (  1 0. 1 FN ) 

WRITE!  10.4419)  PLANE 
4419  FGRMA1  (E15  .7) 
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HR  !  TE  (  IU ,4419 )  F 

HR  1  TE (  1U.442C)  ( C  0 (  I )  .  1  =  1,6) 

4420  FORMAT  t  6  El  5 . 7 > 

C0(b)=CP<6)*180. 

CALL  RELEAS  (5) 

DO  2025  I  =  1  «  NOP 

pm  =p  111*0.72 

Q(  1  )  =  QU>*0.72 
R  (  I  )  '  R ( I  >*0.72 
2025  CONTINUE 
X  5  U  M  =  0 . 

YSUM=0. 

Z  SUM  =  0 . 

DO  20  1  =  1, NGP 
XSUM  =  X  SUM*P ( I ) 

YSUM* YSUV+Q ( 1 1 
ZSUM=ZS UM  +  R  (  1  1 

20  CONTINUE 

00  21  1=1, NOP 

P(  I  )  =  P( I  ) -( XSUM/NQP  ) 

0(  I  ) =U< I  )-( YSUM/NOP 1 
R  (  I  >  =K  (  I  I - < Z SUM/ NO P ) 

21  CONTINUE 

95  I A  =  1 

510  I B= 1 

CALL  OAT A  < C 0 , P , Q , R . F , NOP ,NP AR ,HC ON ,UV LE  , UD  , 
1  VO.  ISP1 ,S0, AM .SPAT ,KM ,UBNO .VBNO, ANOSE) 
TYPE  735.ISP1 
735  FORMAT (  •  ISP1*  •  ,  15) 

WR I T E 1 1 U - 7 1 2 )  ANOSE 
712  FORMAT (E15. 7) 

HR  1  TE  UU  »  71  6 1  KM 
716  FORM  AT ( I  5  ) 

WRITEIIU.720)  i U  BNO (  I ) , VBNDI  I  >  ,  I = 1 , KM  ) 

720  F0RMAT(2cl5.7» 

U5UM  =  0  . 

V  SUM  =  0 . 

DO  10  I=1,15P1 
USUM=USUM+UD(  I  ) 

10  V  S  UM  =  V  SUM* VD (  I  1 
UCG  =  U5UM /  I  SP  1 
VCG  =  VSUM/I  SP1 
WR1TE(1U,782)  UCG.VCG 
782  FORMAT (2E15.7) 

DO  11  I = 1 , I SP 1 
U0<  1 >  =U0 (  I l-UCG 
11  V  D  (  1  )  =  V  l)  (  !  l-VCG 

U2  0=  0 . 

UU2=0. 

U  3  0  =  0  . 


I 

I 

U1 1*0.  - 

003=0. 

U  1  2  --  0  . 

U2  1=0. 

do  so  :  - 1 .  i  spi 

U20=U2Q*UD ( [ ) »-2 . 

U02  =  UG2 ♦  V  0  (  I  I  *»2  . 

U12-U12  +  U0(  i)*vDII)*V0(  i) 

U21-U21+U0I  1  )*Uj(  I  ) *  V  D ( 1  ) 

U i  1  =  L1 1 1  *U0  I  1  *  °  VD  I  1  ! 

U03=U03+VD(  1  i  “>3 . 

20  U30  =  U30*0l)(  ;  1**3. 

U  3  2  U  0  2  /  1  $  P  1 
02  C  =  U  20  /  ISP1 
U3C=u30/ISol 
011=011/  1  S  0 1 
02 1 =U2 1 / ! SPI 

U12=U12/15P1  !- 

U03-U03/ I SP 1 

RM  1  (  I  A  ,  1  6  )  =  02  0  +  U02 

RG-SQRTIRM1 (  I  A, IB)  ) 

RR2I  I  A,  IB)  =  (U20-U0  2  )es>2.*4.*Ull*un 
RH2(  IAvlB)=RM2( 

fi^BI  !A.m  =  (U30-3.*U12)**2.M3.*U21-U03)**2. 

RM3(  IA.  IB)  =  (RM3(  IA.I8)**(4./6.))/(RG**4.) 

RM4  1  IA. IB) = (U30*U12)**2.*<U21*U03>»*2. 
RM4(lA.IB)*(RM4(lA,lB)»*I4./6.))/IRG**4.; 

RM5  (  1  A.  1  B)  *  <U30-3.*U12)  MU30  +  U12)  «*  (  (U30+U12  >  **2. 

1  -3.*(J2i+'JC3)®*2.)*(  3.*U21-U03)»(U2H-U0  3>* 

1  ( 3.  »( J30+U12 .-IU21+U03)* "2 . ) 

R  >1 5  <  I  A.  1  B)  =  (  RM5  (  IA,1B)**I4./12.))/IRG**4.) 

RMfc  l  i  A.  I  B  >  =  (1220-002)  «-(  (  J 3 0*0 12)  ■> ->2  . -I  U 2  !♦  JO  3  >  *»2  . ) 

1  +4  .=01  1-1  U30+U12  )*  (  U2I  *003  ) 

RM6(  IA,  i B )  =  ( R  M6 (  lA,lB)**<4./8.)  >  /  {RG**4 . ) 

0  =  (  2  .  =  'J1 1  )  /  (  U02-U20) 

FBI  I  A,  I B  >  =  A  T  A  N  <  U ) 

FBI  IA,  1 B  )  =  90  .  ®  E  B I  lAf  J8I/3.14 
CU6-EBv!AfIB) 

U20P-0. 

U  0  2  P  =  r' , 

DD  1 26 4  l  -1  .  I  SP  1 

UOD- UO ( 1 )*C0S0(C06>-V0(I )  *S ! NDI 00  6) 

V  0  D  =  U  0  I  I  )®S!ND(C06)  +  Vl»(  1)®C0So(C06) 

U2QP=J20P*UDC’->2  . 
o02P  =  O02P-*V0D  =  :>2  . 

1 264  CONTINUE 

IF  (  U  20  P  .  G  E  .  IJO  2P  )  GO  TO  1276 
! F  (C06.GT.0.)  GO  TO  1278 
C06-C06+90. 

GO  TO  127fe 
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1276  C06-C06-90. 

1276  C0fe'CQ6 

DR  T  -  1 SP 1/150. +  1  . 

1 SP 1 N= 1  NT ( DRT ) 

1  MX  =  Q 

DO  1300  1=1.  ISP1 .15P1N 
1 MX= 1 M  X ♦  1 
1300  CONTINUE 

WRITE!  1U. 13071  IMX 

Du  1562  1  =  1  .  ISP1.1SPIN 

U0(  1  )  =U0  (  1  )  ♦UCG 

VD( ! » =VD( i ) +VCG 

WRITE!  10.1308)  JD (  1).VDII> 

1562  CONTINUE 
1  303  FORMAT < 2E1 5 . 7) 

1  307  FORMAT (  15) 

XI -RM.2 (  !  A, 1 0  ) 

X  2  - RM  3 (  LA, 18) 

X3=RM4(  \  A  «  16) 

X4=RM5 (  I  A, ! B ) 

X5  =  RM6 ( I  A. I B) 

WRITE ( 1U, 1296)  X1.X2.X3.X4.X5 
1298  FORMAT ( 5E1 5 . 7  1 

WRITE  (  III,  1098)  C06.RG 
1098  FORMAT  I 2E 15 . 7  ) 

TYPE  1 198. XI .X2.X3.X4  .X5 
1198  FORMAT ( 51 2X. F7 .4) ) 

TYPE  1398, COB. RG 
1  398  FORMAT  <  2X.F5.1 .4X.F6.4) 

864  CALL  RE  LEAS  l  IU ) 

STOP 

END 


subroutine  data  I  c .P.O.R.F.NCP.NPAR.MCON .uvle.uo, 

1  VO .  I  SP1 . SO  .  AM , SP AT  ,KM  ,UBND .VBND , ANOSE ) 

DIME  NS  ION 

C(6)  ,P(80i  • MCON ( 80,7) ,UBNC(60 ) ,VBND( 60)  , 

1  UD  (  3000) *V0 ( 3000 ) » 1 PN 1 50 ) .0(80) ,R( 8C  ) 

COMMON  KMM .  1UB( 60 )  .  1 VB( 60 ) 

CALL  PROTIC.P.O.R.FjNOP.NPAR.MCON.UBND, 

1  VBND.  KM) 

DO  7  1  =1  .KM 

U6= ( UBNDI 1+1 )-JBNO< 1 ) ) **2. 

V  6  =  (  V  B  N  D  ( 1  +  1 )  -  V  B  N  0  ( 1 ) ) *  *  2 . 

DST1  =  SORT <  Ud*VB) 

PN^-OSTl/UVLE 
IPN  (  1  )  =  1  NT ( PN ) 

PO  =  PN~  1PN (  1  ) 

IF (PO.LT .0.5)  GO  TO  1 


*a  r* 


i 


1 PN< 1 1 = 1PN ( I ) *1 
lFtlPNtl  1  .  NE  .0 )  GU  TG  2 
I  PM  i  >  =  i 

2  DU  = lUBNDt 1 ♦ 1 ) -UBND ( 1 ) ) /] PN (  1  ) 
DVMVbNDt  1+i  )-V6NU(I)  1/iPNI  1  ) 

IF  (  !  . EQ.l)  GC  TQ  4 

I  S°  =  0 

00  3  11=1.1-1 

3  1  SP=  I  SP-rlPNdl  1 
GO  TO  5 

4  ISP  =  0 

5  1PN 1  =  JPNI  I  ) 

00  b  1  I  =1  .  1  PM 
1  A=i  S  P  ♦  I  1 

Ul)(  I  A  1  =  U8ND  (  1)  ♦  <  !  I  -1  1  *0U 
VD  M  A ) =V8ND (I) ♦ (  1 1 -1 > aD V 
CONTINUE 
!  5P  1  =  0 
DU  8  1  =1.KM 

I  5  P 1  =  1  SP 1  ♦  1  PN (  1  ) 
s  continue 

ALNMX  =  -100  . 

ALNMN=100  . 

V2.PAX  =-100  . 

00  1752  1  =  1  .  KM 

IF  (  A  LNM  X . GE .  'JBNO  (  I  )  )  GO  TO  1755 
ALNM  X  =  UBNO (  I  ) 

1755  IF  < ALNMN .LE .UBND! Ill  GO  TO  1752 
ALNMN=UBND< I ) 

1752  CONTINUE 

AL 1  =  ALNMX-ALNMN 
ACNMX  =  -100  . 

ALNMN  =  100  . 

DO  1758  1  =1  , KM 

IF  { ALNMX  ,GE . VBNOt 1) )  GO  TO  1759 
4LNMX=V8NDI I ) 

1  759  IF  I ALNMN  .IE .VBNDl I  ))  GO  TO  175B 
Al_NMN  =  VBND(  I  ) 

1758  CONTINUE 

A C  2  = AlNMX-ALNMN 
TYPE  1751.AL1.AL2 
1751  FORMAT! 2F6. 4) 

A  L  =  (AL1+AL21/2. 

AM  =  0  . 

DO  2147  1=1.1 SP1 
CALL  GAUSSISD.AM.VZ  .SPAT) 

UD (  I  I  =UD ( I  I+VZ 
IF  ( VZMZX.GE  .VZ»  GO  TO  2147 
VZMAX=VZ 
2147  CONTINUE 
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00  2159  1  =  1  .  I  SP  1 

CALL  GAJSSf SO. AM. V2, SPAT ) 
V0<  I  l  *  VO ( I  1+VZ 
IF  (VZMAX.GE ,VZ>  GO  TO  2159 
VZMAX=VZ 
2159  COM  1NUE 

ANCSE=VZVAX*100./AL 
TYPE  1767, ANOSE 
TYPE  1767.VZMAX 
1  7 1> 7  FORMAT  (  F15. 7) 

RETURN 

END 


SUBROUT  1NE  PRDT<  0,P,0,R,F , NOP, N PAR, MCON.UBND. 

1  VBND.KM) 

01  PENSION  0(61, P(BO)  . H  CON  (60, 7),  IJ  (125), VII  25)  , 

1  TA(500).TA1(500).U1NT1(20),V1NT1(20).UINT2(20). 

2  VlNT?(20)tL(20)«LL!2Oi,MKR(2O)>0I!(20)«0I(20), 

3  MKRKC (201  . MX R I  1 ( 20 >  , U6N D<  fcO )  . VBND ( 60 1 , T Au( 6  0 1 , 

9  M C ON  1 ( 12  5*9) •  JK  A  ( 10) tQ( 80)  >R (80) 

D1=C0S0I0I9) ) 

02  =S 1  NO ( 0 ( 9  1  ) 

03s  COSO ( 0 ( 5  1  1 
09  =S I ND( D ( 5) > 

05=  C OS 0 <  0  (  6  )  ) 

06  =  S  I  NO  (  0  (  6  )  ) 

D1  1*05*03 

D12=-05*09*D1+D6*02 

01 3*05*09*02+06*01 

021=09 

022=03*01 

D23=-03*D2 

D31 =-0fa*D3 

D32=D6*09*D1+D5*D2 

G33=-06*D9*02+05*01 

DO  91  1*1, NOP 

PX*P(  l  1*011+0(1  ) =01 2  +  R  <  I)  *013  +  0(1) 

Q  Y  =  P  ( I  1*021  +  011  )*D22+R(  I ) *02  3 ♦ D<  2 ) 

RZ*P<  1 1*031+0  m  *032 +R(  n*033  +  DI3) 

U( 1 )=F*PX/(F+OY) 

V(  1 1 =  F*R  Z / ( F  +  OY ) 

91  CONTINUE 
UMAX*  U( 1 1 
J  X  =  1 

DO  1  1*2. NOP 

IF (UMAX. GE  ,U(  l  1  1  GO  TO  1 


00 


U  N  A  X  =  U  (  1  ) 

J  K  =  I 

1  CONTINUE 
U8N0 ( 1 ) =  U< JK  I 
VBNUI 1 ) =  V< J<  ) 

TAD(l)  =  3H4l5 
K  1=  2 

Ml  =  NOP 
JKD  =  0 

NOP  1 =NQP+  1 
DO  PO  1=1 .NOP 
DO  oO  M  =1  ♦  7 

60  MC0N1  (  I  .  I!) =w  CON  C  !  .  ill 
DO  61  I =  NOP  1  .125 

DO  61  11=1.9 

61  *ICOM  (  1  •  1  1  >  =0 
DO  62  I = 1 .NOP 
DO  62  11=8.9 

62  PCCNl ( I . 1 I) =0 

2  JN1=0 

DO  10  I = 1 • N 1 

IF ( I .EQ.JK  .GR. 1 .EU .JKD)  GO  TC  9 
DO  900  J  =  1  »  9 

1 F ( MC0N1 < JK  .  J)  .Ed  .  1)  GO  TO  401 

400  CONTINUE 
GO  TO  9 

401  1 F  lu(l)  .E3  .U< JK 1 )  GO  TO  3 

xx  =  <  vm-vi  jk i » / c u ( i ) -u(  jk)  ) 

T  A  (  I  1  =A  T  AN  I XX ) 

IF  t J<  I )  .IT ,U( JK)  )  GO  TO  7 
1  F  1  V (  I  )  .IT  .V { JK )  1  GO  TO  8 
GO  TO  10 

3  I  F ( V(  I  ) - V  < JK  )  1  4.5,6 

4  T All ) =4 .7123 
GO  TO  10 

5  T A I  11=20.0 
JN1  =  JN1 +1 

J  K  A ( JN 1  1=1 
GO  TO  10 

6  TA( I )=i .5707 
GO  TC  10 

7  T  A  <  1  )=TAU  J  +  3.1415 
GO  TO  10 

6  TAI  1  1  =T  A  ( I  1+6.2831 

GO  TC  10 

9  TAI  1  ) =20  .0 

10  CONTINUE 

CO  13  1=1 . N 1 

T  A  1 (  I  1  =  T  A  I  I  1-TA0IK1-11+3.1415 
IFITAIIII.LT. 0.0)  GO  TO  11 


1 F ( T A  1 ( 1 ).GE.6.2831>  GQ  TO  12 
GO  TO  13 

11  TA1I  1)=TA1  (  1)*6  .2331 
GO  TO  13 

12  TA1(  I  l-TAl  (  I  )-6  .2631 

13  CONTINUE 
TAM=TA1  ( 1  ) 

J  N-  1 

DO  14  1-2. Nl 

If ( T AH  .  LE  .  TA 1  (  I  )  )  GO  TO  14 
T  A  P  =  T  A 1  (  1  ) 

jn=  i 

14  CONTINUE 

TAO(  Kl )  -  TA(  JN  1 

N2  =  N 1 - 1 
KCLM=0 

DO  26  1 =1 ,N2 
KCL  =  0 

1F( 1 .EO.JK.uR. 1 .  E  Q  .  J  N  )  GO  TO  26 
00  24  11=2, Nl 

1  F ( l I  .EO. JK.OR.I ! .EQ.JNl  GC  TU  24 
1  F (  1  1  .LE  .  1)  GO  TO  24 
DO  500  J  =  1 » 9 

If  (MC0N1  (  I  ,J  )  .EU.l  II  GO  TO  600 
500  CONTINUE 
GO  TO  24 

600  Al=Utll-U(in 
A2  =  J  I JN) -U( JK  ) 

b i - v <  n-vc  m 

B2  =  V ( JN) -V  < JK ) 

A,X  =  A1»32-A2*B1 
ABS AX  =  AB  S ( AX  ) 

IFIABSAX.LE.l .0E-6I  GO  TO  24 

A=<  <U<  JN)-U<  1  1))<-B2-A2*(  V  UN) -VIII  )) )  /  AX 

B=I  AIM V ( JN)-V( 1 1 )) -(U( JN) -U( 1 1)) *B1) / AX 

AL1  =  1  .00000001 

AL2  =  0  .99999999 

EPSl=-0. 00000001 

EPS2=0 .00000001 

1 F< A.GT .All .OR.B .GT .All )  GO  TO  24 
If(A.LT.EPSl.OR.B.LT.EPSl)  GO  TO  24 
IF ( A.LE .EPS2 . ANO .A  .GE .EPS1 )  GO  TU  17 
I F ( A.LE.A11.AND.A.GE.AL2)  GO  TO  18 
GO  TO  19 

17  IF (B.GE.EPS1 .AND .8.LE.EPS2)  GO  TU  24 
1F( B.GE.AL2. AND. B.LE. ALII  GO  TO  24 
KCL=KCL* 1 
L ( KCL ) *  1 

U1NT1  ( KCL)=U(  II) 

VINT1(KCU  =  V(  I  I  ) 


GO  TO  23 

18  1 F ( B  .GE . EP SI . AND .8 .LE .EP S2 )  GU  TO  24 
1FI8.GE.A12.AND.B.LE.AL1 )  GO  TC  24 
KCL=KCL+1 

HKCL  )=2 
UINT1(KCL)=J(  I) 

VI  NT  1 (KCL )  =  V  (  I  ) 

GO  TO  23 

19  IHB.GE.EPS1.ANO.B.LE  .EPS2)  CO  TC  20 
I F { B .GE. AL2.AND .B .LE  .AL  1  )  GO  TO  21 
KCl=KCl+l 

L  (  KCL  )  * 3 

Ul NT1IKCL  >=A*J(  I !♦ { 1 .-A) *U<  1 1) 
VINT1(KCL»=A<*V(  I)  +  (1.-A»*V<  III 
GO  TO  23 

20  KCL=KCL*1 
L ( KCL ) =4 

U 1NT  1  (KCL ) =U( JN  ) 

V 1NT  1  ( KCL ) =V<  JN  1 
GO  TO  23 

21  KCl  =  KCL«-l 
L(KCL)=5 

UINT 11KCL) =U( JK  ) 

VINT1(KCL)=VIJK) 

23  MKR ( KCL  )  =  I  1 

A  A=  U  (  JK)  -UINTK  KCL) 

BB=  V  < JK)-VINTI(KCL) 

D1  1 (KCL) =  5QRT ( A A* *2 .♦88’ *2  .) 

24  CONTINUE 

IF( KCL.EQ.O)  GO  TO  26 
kclm  =  kclmm 

MKR  KC (KCLM ) =  1 
01  (KCLM) =01  1111 
Jl  =  l 

IF ( KCL  .EG  .1  )  GO  TC  5C 
00  25  III =2. KCL 

I F ( 01 (KCLM )  .LE  .0  I  1 (  I  I  I  ) )  GO  TO  25 
0  I ( KCLM ) -0  I  I  (  I  I  I  ) 

J  1  =  I  I  1 

25  CONTINUE 

50  MKR l  I ( KC  LM  )  =MKR ( J 1  ) 

LL ( KCLM) =L ( J  1  ) 

UINT2I KC  LM  )  =  U  1NT 1 ( J1  ) 
VINT2IKCLMJ-VINT1 ( J1  ) 

26  CONTINUE 

I F ( KCLM.EO.O)  GO  TO  34 

DM  N -  0  l  (  1  ) 

LPA  =  1 

1 F ( KCLM  .EQ.ll  GO  TO  28 
DO  27  1=2. KCLM 


1  F  (OWN  .LE  .0  l  (  1  )  )  GO  T(J  27 
DMN1 D  l  (  I  ) 

LPA-  I 

27  C  ON  T l NUE 

28  UB  NO  < Kl ) =01 NT2( LPA) 

VBND(Kl) =VI\T2(LPA) 

L I =MKRKC (LPA) 

LlI=f*KRIl  (LPA> 

DO  100  J*1.9 

I F (MC0N1  ( L  I  .  J  )  ,EQ  .L  l  tl  GO  T  0  200 
GO  TO  100 
20C  J  1  1  - J 

GO  TO  300 

100  CONTINUE 

3  OC  00  101  J-1,9 

I  F(  MCCM  (  L  I  1  ,  J)  .EQ.L  1 )  GO  Ttl  201 

GO  TO  101 

201  J  1 2  =  J 

GO  TO  301 

101  CONTINUE 

301  00  1020*1.9 

I F ( MCON 1 ( JK , J ) . E 0 . ON )  GO  TO  202 

GO  TO  102 

202  J 1 3  —  J 

GO  TO  302 

102  CONTINUE 

302  DO  103  J*l»9 

I  F ( MC0N1 ( JN  .  J I . EQ . JK )  GO  TO  203 
GO  TO  103 

203  J 1  A  =  J 

GO  TO  303 

103  CONTINUE 

303  I F( L L ( LPA ) .EO  .1  )  GO  TO  29 
I F ILL ( LP A)  .EO  .2  1  GO  TO  30 
IF(LL(LPA).E0.3)  GO  TO  31 
1 F ( L  L ( LP  A  5  .  tQ  .  4 )  GO  TO  32 
MC0N1 ( L 1 . J  Ill  * JK 

NC0N1  (LI  1 , J12  )«JK 
MCON  1  ( JK  ,8) =L  I 
NC0N1 ( JK, 95  *L I  I 
GO  TO  2 

29  MCCM  (  JK,  Jl  3)  =L  1  1 
NCCN1  ( JN.J  1  4  )  =  L I  I 
MCCN1  (LI  1  ,85  =  JK 
KC0N1  !L I  1  ,9  )  = JN 

J  K  D  =  J  K 
J  K  =•  L  II 
GO  TO  33 

30  VC0N1(JK,JI3)=L1 
MC0N1  ( JN . J I  4 ) =L  I 


V  c  DN 1  ( L I .8 ) * JR 
MC  ON  1 ( Ll  •  9  )  s  JN 
J  kO^jk 

jkh! 

GO  TO  33 
N  1 1  N  1  ♦  1 

MCCN1 (JK.JI3) *N1 
MCCN1 ( JN  ,J 1 4  1  *N1 
MCON  1  <  L  I  .  J  II  )  =  N1 
MC0N1 { Li  l . J 1 2  )  *N1 
MCON l ( N1  . 1  Is  L  I 
MCCN1 «  N  1 . 2 ) =L  l  1 
MCON1 <  N1 .3) *JK 
MCON 1 (  N1 ,4)= JN 
MCON1 { N 1 . 5  )  -0 
MCON 1 ( N1 %6  )  *0 
MC0N1 (N  1 , 7  )  *0 
U ( N 1 ) -UB  ND ( K  1 ) 

V ( N 1 ) =  V  6  ND ( K 1 ) 

J«D  =  JK 
JK  =  N  1 
GO  TO  33 

MC0N1 ( L  l  . J 1 1 ) 1 JN 
MCGNHL1  I  ,J12)*JN 
MCON  1  <  JN  .8  >  *U 
MC0N1 ( JN.9  )  *L  II 
JKD= JK 
JK  =  JN 
K  1  =  K  1  ♦  1 

^(U?JN»  .EO.UBNOd  )  .ANO.VUN)  .EQ.VBNDim 

1  GO  TO  35 
UBND ( Kl )  =U( JN) 

VBNOIKl  )  =V( JN) 

JKD- JK 
JK  =  JN 
K  1  =  K  1  ♦  1 

GO  TO  2 

KM  =  K 1-1 

UBNut  KMfl ) =  U  BN  0 ( 1 ) 

VBNO ( KM+1 ) = VBNO< 1 ) 

RETURN 

END 


SUBROUTINE  GAUSS < 50.AM  ,  VZ  .SRAT  ) 
A  =  0 . 
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Y=R AN ( SPAT  I 
50  A  =  A  ♦  Y 

V Z  =  ( A-fe  .  ) *SD+AM 

RE  T JRN 

END 


DIMENSION 

X11(19,19)sX12<19,19).X13(19,19),X14(19,19), 

1  X15( 19.19)  iRKl ( 19, 19) , E  B 1 ( 19. 191 
DIMENSION 

X21(19, 19). X22! 19,19)  .  X  2  3  ( 19,19)  » X  24  (  1  9  # 1 9  )  , 

1  X25( 19.19) . RK  2 ( 1 9 , 1 9  ) ,  EB 2  ( 19. 19) 

DIMENSION  X1(11,11),X2(11,11),X3<11,11), 

1  X4<  11. 1 1 ) . X  5 (  1  1  ,  11  ) 

DIMENSION  CO(fe) 

DIMENSION  UDI220),VO<220).UtiNG<60),YBND(feO> 

TYPE  100 

100  FORMAT! '  ENTER  3  C  TAPE  NC.  ') 

RE AD i5. 120)  lU 
120  FORMAT ( I  1  ) 

1  U  = I U  +  8 

I  F  !  IJ.NE.17)  GO  TO  140 
TYPE  150 

150  FORMAT!'  LGAOING  OF  MQNENT  TABLES  ') 

TYPE  170 

170  FORMAT!'  ENTER  DEC  TAPE  NO.  ') 

READ! 5, 190)  IU 
1 90  FORMAT (11) 

I U* 1 U ♦8 

CALL  I  FI  LE (  lu,  'PHAN '  ) 

DG  610  I  -1  ,  19 
DO  fell  11=1.19 

READ! 1U.250)  AZ.RQ.X11(I.II),X12(I,I1).X13(I.II). 
1  X14 (  I  . 1  I  ) . X 1 5 (  I  ,  I  I  )  ,RK1 (  I  ,  I  1  )  , EB1 (  I  .  II) 

611  CONTINUE 
610  CONTINUE 

CALL  RELE  AS  (  IU  ) 

CALL  IF  I  LE  (  1 U. 'M  Ift  '  > 

DO  310  1=1.19 
DO  31  1  11  =  1.19 

RE  AD  (IU. 250)  AZ.R0.X21(I,m,X22II,m,X23(I,JI), 
1  X24!  I, I  I  )  ,X25(  1,1  I  )  »R*  2  (  I  .II  )  ,EB2(  I  .11  ) 

311  CONTINUE 
310  CCNTINUE 

CALL  RELEAS  ( 1 J) 

25j  FORMAT (2F6.1.7E15. 7) 
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TYPE  390 

390  FORMAT! '  END  OF  LOADING  ') 

R  E  AD ( 5 , 2  50  )  ENO 
KO  IF  (  10.60.8  )  GO  TO  <,20 
TYPE  430 

430  FORMAT!'  ENTER  FILE  NAME  AS  »**  ',/) 
READI5.450)  IFN 
450  FORMAT ( A3 ) 

C  AIL  IF  ILE!  IU.IFN) 

READ! 1U.480)  PLANE 
480  FORMAT! E15. 7) 

READ!  IU.480)  F 

RE  AD (  IU,  500  I  !C0! I  ), 1  =  1 ,6) 

500  F0RMAT(6E15.7> 

READ!  IU.480  )  ANOSE 
READ!  IU, 530  )  KM 
530  FORMAT (15) 

READ! IU. 550)  ( UBND (  I  ) , VBND (  I ) ,  I =1 , K M ) 

350  FORMAT (2E15. 7) 

READ! 1U.550)  UCG.VCG 
READ!  IU  ,  530 )  ISP1 

READ ( IU, 550  )  (UD  <  I  )»VO(l).I  =  l.ISpl) 

READ! I U . 600 )  Y1 , Y2  ,Y3 , Y4, Y5 
600  F0RMAT(5E15.7) 

READ!  IU .  550 )  C06,RG 
CALL  RELEAS  <  IU) 

420  DM  1  - 1 000  . 

DM  2  =  1 000  . 

TYPE  635, ANOSE 

635  FORMAT!  •  NOISE*  '  ,F6.1) 

TYPE  637, PLANE 

637  FORMAT!'  AIRCRAFT  CUDE*  ',F3.1> 

CO!  1>=C0(1)/12. 

CO ( 3 ) *C 0 ( 3 ) / 1 2  . 

CD(2)=C0(21/12. 

TYPE  1190. CO!  1)  . COO) 

TYPE  1210. CQ(2) 

TYPE  1230,CO(4),CO(5),CO(6) 

CALL  DI S IUBND.VBNO  ,KM ,UD,VO , ISP1) 

DO  640  1*1,19 
DO  650  11*1,19 

GR1=  (  Yl-Xl  1  (  1  . 1  1  )  )  4*2. •»  (  Y2-X12!  I  ,1  1  )  )  «*2. 
1  ♦ ( Y3-K13 (  I  .  I  I  )  )  **2 . ♦< Y4-X14 (  I  , 1  I  )  ) »*2. 

1  ♦ t  Y5-X15!  1 .  1  I )  ) **2 . 

IF  (DM1  .LE . DR  1 )  GO  TO  690 
DM  1  =  DR  1 
I  A  1  =  I 
I  R  1  *  I  I 

6  90  DR2  -  (  Y1-X21  (  I  .  I  I  >  )  **2  .♦(  Y2-X  22!  I  » I  I  )  ) **2  . 
1  ♦ ! Y3-X23!  I  .  I  I  )  ) **2 .♦( Y4-X24 (  I ,  II)  ) **2  . 


1  ♦  <  Y  5-X  2  5  C  1  ,  I  I)  )  **2. 

IF  (DM2.LE.DR2)  GO  TO  650 
I  A 2=  1 
IR2=1  1 
0*2  =  DR2 
<,50  CUNT  1  MUE 
640  CONTINUE 

IF  (DM2. Lc. 0*1)  GO  TC  820 
CER=DM2/ (DM2+DM1 ) 

P  l  =  1 . 

CQ6=C06-E31(  I  A 1 . 1 R 1 > 

C02  =  BK1 (  I  Al. I R 1  l/RG 
CALL  1  N  T  P  (  IA1  •  I  R 1 1  XU  .  X 1 1 
CALL  INTP( IA1 .IR1.X12.X2) 

CALL  INTPI  I  A  1  .  I  R 1 .X13, X  3 ) 

CALL  |NTP(  IA1 .IR1.X14.X4) 

CALL  INTPI  IA1.U1  .X15.X5) 

GO  TO  915 

620  CER=DM1/DM1+DM2 
PL  =  2 . 

C06  =  C06-EB2  I  I A2, 1 R2) 

CC2  =  RK2(  IA2 , IR21/RG 
CALL  INTPI 1A2.1R2.X21, XI) 

CALL  INTPI  I A2  .IR2.X22.X2) 

CALL  INTPI  1A2.IR2,X23,X3) 

CALL  INTPI I A2.IR2.X24, X4) 

CALL  INTPI I A2 . 1R2 , X25,X5 ) 

915  CER'CER  *100 . 

DM  =  1000 . 

DO  1010  1*1.11 
DG  1020  11*1.11 

D  =  < Y  1-X1 1  I  ,  I  II ) **2  .♦( Y2-X2 I  1  .  I  H ) **2  . 

1  +IY3-X3I  I.  1  1)  1**2. ♦( Y4-X4 11,11)1**2. 
1  -.IY5-X5I  1  .  1  1)  1  *  *2  . 

IF  ( DM  .  LE  .0  >  GO  TO  1020 
DM  =  0 
I  AD  =  I 
1RD*  I  I 

1020  CONTINUE 
1 C 1 0  CONTINUE 

C01=UCG* (F+CG2) / F 
C03  =  VCG* I F+CC2) /F 
TYPE  1122 

,.122  FORMAT  I  '  ESTIMATED  VALUES  ’> 

IF  (PL.NE.1.1  GO  TO  11 3J 
C05  =  5.*IIA1 -1.1*1  1  AD-6.  ) 

C  04  =  5  .*(  1R1-  1  . )  ♦  I  IRD-b  . ) 

TYPE  1355.IA1.1R1.IAD.1RD 
1355  FORMATI4110) 

TYPE  1 1 52 ,CER 
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1152  FORMAT ( 1  AIRCRAFT  IDENTI F! CO  TC  BE 
PHANTOM, CERTAINTY- ' , 

i  fs.ii 

GO  TC  1160 

1130  C05*5 .* (1 A2-1 . ) ♦ ( I  AO-6. ) 

034  =  5  ( 1 R2-1 .  ) ♦< I RD-6  .  1 

TYPE  1355. .A2.1R2, 1A0.1R0 
TYPE  1162, CER 

1162  FORMAT ( '  AIRCRAFT  IDENTIFIED  TO  BE 
MIRAGE .CERTAINTY* ' . 

1  F5.1) 

1 1  (>C  CGi=CUl/L2. 

C Q 3 * C 03/ 1 2  . 

C02*C02/ 12  . 

TYPE  1190.C01.C03 

1190  FORMAT!'  X*  ' , F6 . 1  •  '  Z*  '.F6.ll 
TYPE  1210. C02 

1210  FORMAT ( '  Y*  '.F10.ll 
TYPE  123G,C34.C05,C0fe 

1  2  30  FORMAT!  '  ROLL  ANG  .  =  ',F5.1,'  AZiMJTH  A NG  .  =  '. 
1  F  5 . 1  « '  ELEVATION  ANG  .  =  '.F5.ll 
STOP 
END 


SUBROUTINE  INTP( I  A  .  I  R  »  X  X  ,  T  1 

DIMENSION  XX( 19.191  ,T(  11,111 

A  5  *  X  X  (  I  A  .  IR  1 

IF  (1A.E0.11  GO  T  G  30 

A  2  s  X  X (  IA-1  ,  IR  I 

IF  I  IR.NE  .1  1  GO  TO  40 

Al* 1000 . 

A 4  --10  0C  . 

A7* 1000  . 

A  6  =  X  X  (  1  A  ,  I  R  ♦  1  1 
A  3  =  X  X  (  I  A  -  1  ,  I  R  ♦  1  1 
GO  TC  50 

40  A1  =  XX(  1  A-l  ,1  P.-l) 

A4  *  X  X (  1A.1R-1) 

IF  (1R.NE.19)  GO  TC  60 
A3* 1CJ0 . 

Afc*  1000  . 

A9 *  1000  . 

GO  TC  70 

6C  A3=XX( 1 A-i , IR+1) 

A  6  =  X  X  (  1  A  ,  I  R  +  1  1 
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IF  (1A.NE.19)  GO  TC  80 
A  7  -  10  JO  . 

A8- 1030  . 

A9  =  1000  . 

GO  TO  530 

80  A7  =  XX ( 1 A»1 , I R  — 1 1 

A 8  =  X X ( I  A* 1  •  l R ) 

A9  1  X  X  (  1  A«-1  .  1  P  *1  ) 

GO  TC  500 

50  IF  (  IA.NE.19)  GO  TC  100 
A8  =  1 0 00  • 

A 9= 1 000  . 

GO  TO  500 

ICQ  A  8- X  X  (  I A+ 1  .  I  R  ) 

A9  *  X  X  (  1  A+  1 1  1  R  -*■  1  ) 

GO  TO  500 

70  I F  (  !A  .NE  .  19  )  GO  TC  120 
A7=  1000  . 

A8= 1030  . 

GO  TO  500 

120  A7  =  X  X ( 1 A+l . 1R-1 > 

A6  =  XX ( I  A* 1  .  1R  ) 

GO  TO  500 
30  A 1  =  1 000  . 

A2= 1000  . 

A3  =  1000  . 

A8  r  X  X ( 1  A ♦ 1  .  1R) 

IF  ( 1R.NE .1  )  GO  TU  200 
A4=  1000  . 

A7* 1000  . 

A  6  -  X  X  {  1A.IRM) 

A9  =  XX I1A+1.  1R*1  ) 

GO  TO  530 

200  IF  (Ift.NE  19)  GO  TO  300 
A6= 1000  « 

A9=  1000  . 

A  4  =  X  X  !  A,  1R-1  ) 

A 7 s  X X I 1 A*1  ,  1R“  1  ) 


GO  TO  500 
A4 -XX ( 1  A  ,  1 R-  1  > 

A6  =  XX (  1  A  .  1R+  1 1 
A  7  =  X  X  (  l  A  ♦  1 .  IR-l  > 

A9 1  X  X  (  I  A+l  ,  IR+1  )  ,  ..  .a  n 

CALL  EQPT( A1.A4.A7 .T (1.1  ),T(2,l).T(j.l).Tt4.1), 

r(5,l)«T(fc.l)iT(7.1).T(8,l).T(9»l).T(10.1).T(ll.l)) 


300 


5C0 


EGPT(A3*Afe.A9.T(l.ll).T(2.11).T(3.11)%T(4.11)» 
T(5.11).T(g.U>.H7.11)  ,T(8,11)  .T(9,ll)  .T(10,ll). 


1  TUI, ID) 

CALL 

EGPTIA2,AS,A6,T(1,6).T(2,6).T(3,6),T(4,6),T(5,6), 

1  T  (  6 . 6  )  ,T(7,t).T(8,6)  *  T  (  9 , 6  )  ,  T  (  1  0  *  6  )  , 

1  T ( 11.6  )  ) 

CALL  EQPT (T( 1  ,1  )  ,T ( 1,6 ) ,T< 1 . 11)  , 

1  T  (  1,1)  ,T(  1 , 2 ) , T ( 1,3) ,T(1.4),T( 1,5)  , 

1  T(1,6),T(1,7),T(1,8).T(1,9),T(1,10).T(1,11)) 
CALL  EQPT< T< 2 . 1 > ,T (2,6) ,T(2, 11 ), 

1  T (2.1)  *T (2,2)  ,T( 2.3) , T ( 2 , 4 ) .  T  (  2 . 5  ) , T  t  2 , 6 )  , 

1  T(2,7),T(2,8).T(2,9).T(2,10),T<2,11)) 

CALL  EQPT  !T( 3.1) ,T(3,6)  ,  T ( 3  . 1 1 )  . 

1  T(3.1>.T(3,2).T(3,3),T(3,4»,T(3,5).T(3,6>, 

1  T(3.7)»T(3,8),T(3,9)»T(3,1C),T(3,111) 

CALL  EQPT(T(4, 1) ,T<4,6) ,T(4,11)  , 

1  T(4,1),T(4,2)»T(4,3)»T(4»4),T(4,5),7(4,6), 

1  T(4,7),T(4,8),T(4,9)*T(4,10).T(4,11)) 

CALL  EQPT(T(5.1).T(5.6) . T  (  5 , 1 1 )  , 

1  T(5.1),T(5,2).T(5,3)«T(5.4)»T(5»5),T(5,6), 

1  T(5,7),T(5,8).T(5,9) . T ( 5 , 1 0 ) , T (5,11) ) 

CALL  E0PT(T(6«1),T(6,6)«T(6,11), 

1  T(6.1).T(6,2),T(6,3)*T(6,4),T(6,5),T(6.6), 

1  T(6«7),T(6»8)»T(6*9),T(6,10>,T(6,11)) 

CALL  EGPT( T< 7. 1) ,T ( 7,6) ,T( 7. 11) . 

1  T  <  7 . 1 )  ,T(7,2).T(7,3) ,T(7.4),T(7,5),T(7,6), 

1  T (7,7),T(7*8),T(7»9) »T (7,10),T  (7,11)  ) 

CALL  EOP T ( T < 8. 1) ,T ( 8  ,6)  ,T<8  , 11)  , 

1  T(8,1),T(8,2),T{8,3),T(8,4),T(6,5),T(6,6), 

1  T(8.7),T(8,8).T(0,9)»7(8,1O),T(8,11)) 

CALL  EOPT(T(9,l) ,7(9,6)  ,T19, 11) , 

1  T(9.1),  7(9, 2). 7(9, 3), 7(9. 4), 7(9, 5), 7(9. 6), 

1  T (9.7) ,7(9,8) ,7(9,9), 7 ( 9 , 1 0 ) . 7 (9,11  ) ) 

CALL  EQPT(7(10,l).T(10,6)rC(10,ll), 

1 

T(10, 11.7(10, 2). T(10, 3), 7(10, 4), 7(10, 5), 7(10,6), 

1  T ( 10,7)  ,7( 10,8) »T( 10,9) ,T(  10,10)  ,7(  10,11)  ) 
CALL  EQP7(7< 1 1,1 ) ,T(11 ,6) ,T(  11  , 1 1  )  , 

) 

Till, 1;, 7(11, 2), 7(11, 3), 7(11, 41, 7(11, 5), 7(11, 6), 

1  7(11, 7),  7(11, 8),  T(ll,  9),  7(11, 10),  T(ll,  ID) 

RETURN 

END 


SUBROUTINE  EQPT(X1,X2,X3,Y1,Y2,Y3.Y4,Y5.Y6,Y7,Y8 
1  Y9,Y10,Y11) 

01 = ( X2-X1 ) /5 . 
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D2  = (X3-X2) / 5  . 

Y  1  =  X  1 

Y  2  =  Y  1 +  01 
Y3  = Y2  +  D1 
Y4  *  Y  3*0  1 
Y5=Y4+D1 

Y  b  =  X  2 
Y7=Y6*D2 
Y8=Y7*02 
Y9-Y8+D2 
Y10-Y9+Q2 

Y  1  1  =  X  3 

return 

END 


SUBROUT  1 NE  D I S ( U BND , V  BN D . KM , UD . V D . 1 SP1  > 

DIMENSION  UBNDI 60 >  .VBNDI 60) 

DIMENSION  UD1250 1  ,VD(  250) 

DIMENSION  A ( I  000 ) 

UBNDIKM+1)  =  UBND(  1  I 
VBND(KM*1)=VBND( 1  ) 

C  0= 1 0000 . 

COMMON  LP.ISHOW.XMAX.XMIN.YMAX.YMIN.INTENS. 

1  1  SC  ALE 
LP  =  0 
I  SHOW  *  0 
INTENSE 
I  SC AL E  =  0 
KMM-KM>1 

13  DO  1  I *1,KMM 

UBNDI I l*CD*ueND( 1 MlOOO. 

VBNDI  I ) =  CD*V4  NDI  11*1000. 

I  CONTINUE 

DO  65  I  =1  .  1  5P1 

UD ( I >  =  C  D ®UD ( I  1*1000. 

V  D I  1 ) =CQ«VD( 1 1*1000. 

65  CONTINUE 

II  CALL  INTAB  I  A, 1000) 

XM  AX  =  2000  . 

X  M  1  N  =  0  . 

YMA  X  =  2000  . 

Y M IN  =  0  . 

CALL  POINT  U, UBNDI  1  )  .VBNDI  1  )  1 
I  SHOW  =  1 
DO  2  I  *1 ,KM 

CALL  LINE! A.UBNOII  ) , V BN D I  I  I  * UBNDI 1*1)  .VBNDI 1  +  1) ) 


2  CONTINUE 

CALL  PCI  NT ( A, 53  .  .5C  .  ) 

CALL  LINE ( A.  50  ..  50  . ,  1950  .,50. 1 
CALL  LlNE(A,1950.,50.«1950.»  1950  .  ) 

CALL  LINE ( A. 1950. ,  1950. ,50. .  1950  . ) 

CALL  L I NE ( A. 50  . ,  1  950  . ,50  .  ,50  .  ) 

DC;  35  1  =  1  .  1  SP  1 

CALL  SPOT ( A.UU!  I) ,VD< I  >) 

35  CONTINUE 

CALL  DlSPLY(l.A) 

TYPE  3 

3  FORMAT!'  TYPE  MAG.  FACTOR  IN  FORMAT  F4.1 
R  E  A 0  <  5 . 7 )  AB 

7  FQRMATIF4.il 

15  DO  06  1  =  1  ,  1SP1 

UD<  1  )=(U0(  I  1-1000. 1/CO 

v o ( i  )=(vo(  n-iooo.i/co 

fefc  CONTINUE 

00  3252  1=1. KMM 
UBNO ( 1 ) = ( UBND ( 1 1-1000.1 /CD 
V  BN 0 ( I )  =  ( V  BND ( I  1-1000 . 1 /CO 
3252  CONTINUE 
CD=  A6»CD 

IF  (AB.NE.l  ,1  GO  TO  13 

RETURN 

END 
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